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

УДК 531.383:532.516

Математическое и компьютерное моделирование нелинейных волн деформаций в оболочке, содержащей

вязкую жидкость

Ю. А. Блинков*, С. В. Иванов*, Л. И. Могилевич^

* Механико-математический факультет Саратовский государственный университет им. Н.Г.Чернышевского г. Саратов, Россия, 410012 ^ Поволжский филиал Московский государственный университет путей сообщения г. Саратов, Россия, 410790

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

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

1. Введение

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

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

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

Кроме того, знание базиса Грёбнера дает возможность проверить совместность исходных разностных соотношений, определить произвол в решении, посчитав полином Гильберта, и, применяя специальный вид допустимого упорядочения при его построении, получить другое представление первоначальных разностных соотношений.

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

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

Работа поддержана РФФИ (гранты 10-01-00177-а «Задачи упругогидродинамики ребристых и трехслойных пластин и цилиндрических оболочек» и 10-01-00200-а «Моделирование квантовых вычислений и дискретных систем методами компьютерной алгебры»).

2. Постановка задачи гидроупругости и построение математической модели методом возмущений

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

дУ- „дУ- дУ- 1др [д2У- 1дУ- д 2У- У- \ —- + К—- + + ~ = V + +

д

дК дг

д

дх рдг \ дг2 г дг

д х2

, т. дУх дУх 1 дР (д2Ух 1 дУх д2Ух + Уг—--+ Ух—--+ - — = V I +----—+

д

дх р дх

дУ- У- дУх п

1Г + = °

д д х

д 2 д

д х2

На границах с оболочками выполняются условия прилипания жидкости

(1)

У- = -

дШ

У = дV

при

г = КЛ-Ш.

(2)

Здесь Ь — время; Уг,УХ — проекции вектора скорости жидкости на оси цилиндрической системы координат; р — давление; р — плотность; V — кинематический коэффициент вязкости; V — продольное упругое перемещение оболочек по оси х; К — прогиб, положительный к центру кривизны оболочки; К\ — внутренний радиус оболочки.

Записывая уравнения движения элемента цилиндрической оболочки в перемещениях для модели Кирхгофа-Ляве, считаем материал нелинейно-упругим с кубической зависимостью интенсивности напряжений а\ от интенсивности деформаций е1 [4,5]

а1 = Ее 1 — т е(3)

Здесь Е — модуль Юнга; т — константа материала, определяемая из опытов на растяжение или сжатие.

Уравнения динамики физически нелинейной оболочки с учетом (3) записываются в виде [5]

ЕН° д_ 1 — ^2 дх

'ЗЦ_ д х

_ 4 тН° д_ № к) — 3 1 — р2дх

,диШ (ди

— (1 + №) ЖкЫ — я) —

ди

д х

К

ЕН° (коди _4 тЬ° 1

1 —р20\ 12"дх4 — ~К~дх + В?) — 3 1 — р2° К

— р°Ь°

д2 и дг2

( - У

д х

(4)

.. .диШ [ди К . ...

— (1 + ы -вххвУдх — К) — IК

)—(К)

+ р°Ь°

д2 К

3

3

х

3

а

Здесь р° - плотность материала оболочки; — коэффициент Пуассона; К — радиус срединной поверхности оболочки; Ь° — толщина оболочки (Н°/2 = К — К1);

Со = \/Е/[р0(1 — ^о)] — скорость звука в материале оболочки; дп — напряжения со стороны жидкости.

Если снести напряжения на невозмущенную поверхность оболочки (Щ ^ К), то можно считать, что поверхностные напряжения со стороны жидкости определяются формулами

<1х =

рр

/ ЗУ, + &К \ \ дг дх )

Яп =

г=Н

+ 2 дУг

—р + 2р и—-д

(5)

г=Н

Принимая за характерную длину — длину волны , перейдем к безразмерным переменным для исследования уравнений (4)

х

Щ = тти3, и = итщ, г* = у^, х ^

(6)

Положим

^ = е = о(1), ^ = 0(е),

К = оа

К Е

К = 0(е1/2),

К

Е = 0(е), е9 К = 0( ел/2).

(7)

т

где е ^ 1 — малый параметр в задаче (4).

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

£ = х* —сГ, т = еГ

(8)

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

и1 = ию + еиц + ... из = изо + еиз1 + ...

(9)

Для определения правой части уравнения (4) ведем безразмерные переменные и параметры

Уг = Шт — Уг, Ух = Шт — Ух, Г = —, г = — г, х = -, р

К1 К1

ррсо1шг,

Кз1

Р,

К1

ф = К1 = 0(£ 1/2), Л = Шр = 0 (е), ф « 1, Л < 1.

К1

(10)

Подставляя (10) в уравнения гидродинамики (1) и граничные условия (2), представим безразмерные скорость и давление в виде разложения по малому параметру Л

ух = + Л VI + ..., уг = + Л VI + ..., Р = Р0 + ЛР1 + ... (11) С принятой точностью по Л,ф,£ положим К ~ К и окончательно получим

д2ию , 1 2 1 — р20 д4ию 2т /ит\2 , Г 2

+ Дт) -2--д? — жЫ ^ — № + V1 — ^

/ дщЛ V д^ )

д2ию 2

2

1 — (2№ )2

р д и

10

роЪ,о£В,1Со

= 0. (12)

2

х

Легко видеть, что замена

ди

10

д

= ср, г/ = С1 Ь = С2Т

позволяет записать уравнение (12) в виде

дР + д!р_6 2др_а =о

дЬ дг]3 дг]

(13)

(14)

Здесь а = +1 (при ро < 1/2 — неорганические материалы), а = _1 (при ро > 1/2 — живые организмы) и а = 0 (при ро = 1/2 — резина).

Постоянные с, С1, С2 определяются при подстановке (13) в (12) и имеют вид

С2 = а 2

1 _ (2ро)2

р1р

роИое^со Е К2

1

2

2

К) Ро\/1 _Ро

1/3

=

Ро

2т и2^ 1 _ ро + Ро_

1/2

С1.

В случае отсутствия жидкости последнее слагаемое в уравнении (14) исчезает, переходит в МКдВ и имеет точное частное решение в виде кинк-антикинк

р = ±к Ып(к(г] + щ) + 2к3 £).

(15)

Эти решения при £ = 0 можно взять в качестве начальных условий при решении задачи Коши для уравнения (14).

(

3. Компьтерное моделирование

В работах [6-8] развит подход к построению разностных схем, основанный на построении переопределенной системы разностных уравнений, получаемой из аппроксимации интегральных законов сохранения и интегральных соотношений, связывающих искомые функции и их производные. В результате разностная схема определяется как условие совместности для данной системы. Таким образом получается разностная схема метода конечных объемов, автоматически обеспечивающая выполнение интегральных законов сохранения по областям, составленным из базовых конечных объемов.

Запишем уравнение (14) в интегральной форме

/(2р3 _ р^+ рйг) _ = 0 (16)

■ 'п

дП

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

п

У щ&г] = и(г, ^+1) _ и(г, щ),

(17)

Щ+1 К '

J ичп¿г] = и.п(г,щ+1) _и.п(г, щ).

т

n + 1

j j+1 j+2 Рис. 1. Базовой контур для уравнения (16)

n

Используя для интегрирования по времени и по четным производным по г] формулу трапеций, а по нечетным производным по г] формулу среднего значения, и полагая гп+1 — гп = т, щ+1 — щ = И, перепишем соотношения (16), (17) в виде

2и3" + 2и3 п+1 — 2и3"+2 — 2«3п+1) —

п I „. п+1 _ „ п _ „ п+1\\ _

[y-VVj + UVVj UVVj+2 UVVj+2 )) ' 2 +

,n+l „j \ Oh

+ (u^1 - un+1) ■ 2h - a(un+-t + un+1) ■ hr = 0,

h

iuvj+i + uvj) ' 2 = Uj+1 — uj ,

U п ■ 2h — u п — u п uVVj+1 2h — uVj+2 uVj ■

Используем пакет [9] для построения базисов Грёбнера в случае разностных идеалов [10], реализованный в системе компьютерной алгебры Maple [11].

> restart:

> libname:=libname, "/usr/local/lib/LDA"":

> гсг^^БЛ):

> Ь: = [((2*Р(П^)+2*Р(П+1^)-2*Р(П^+2)-2*Р(П+1^ +2)) -

> (гак^^) +uxx(n+l,j)-uxx(n,j+2)-uxx(n+l,j+2))) *tau/2+

> (u(n+l,j+l)-u(n,j+l))*2*h - sigma*(u(n+l,j+l)+u(n,j+l))*h*tau,

> (ux(n,j+l)+ux(n,j))*h/2-(u(n,j+l)-u(n,j)),

> 2*uxx(n,j+l)*h-(ux(n,j+2)-ux(n,j))];

Ь := [1/2 (2Р (п,¿) + 2Р (п + 1,з) — 2Р (п,з + 2) — 2Р (п +1,^ + 2) — ихх (п,з) — ихх (п + 1, з) + ихх (п,з + 2) + ихх (п + 1, з + 2)) т + 2 (и (п + 1,з + 1) — и (п,з + 1)) И — а (и (п + 1,з + 1) + и (п,з + 1)) Ит, 1/2 (их (п,з + 1) + их (п, ])) И — и (п,з + 1) + и (п, у), 2 ихх (п,з + 1) И — их (п, j + 2) + их (п, з)]

Поскольку пакет [9] работает только в случае линейных разностных идеалов, а исходное дифференциальное уравнение (14) нелинейно, заменим нелинейную часть введением дополнительной функции Р. За счет выбора допустимого лексикографическое упорядочение сначала по функциям ихх У их У и У Р, затем по переменным п,з, нелинейная часть не будет входить в лидирующие мономы системы при построении базиса Грёбнера.

> JanetBasis(L, [П^], [uxx,ux,u,F],2);

[[-г и (п, j) + 2ru (п + 1, j + 1) + (4 h3 - 2 Н3ат)и (п + 1,j + 2) -ти (п + 1, j) + 2ru (п, j + 1) + (-4 h3 - 2 h3aт) и (n,j + 2)

- 2т и (п + 1, j + 3) - 2 г и (п, j + 3) + т и (п + 1, j + 4) + т и (п, j + 4)

- 2 h2T F (п + 1, j + 3) + 2 h2T F (п + 1, j + 1) - 2 h2T F (п, j + 3) + 2 h2r F (п, j + 1) hux (п, j + 1) + hux (п, j) - 2 и (п, j + 1) + 2 и (п, j), h2uxx (п, j + 1) - и (п, j + 2)

+ 2 и (п, j + 1) - и (п, j), -Ь2тихх (п + 1, j) - Ь?тихх (п, j) + т и (п + 1,j + 3)

- 2т и (п + 1, j + 2) + (4 h3 - 2 h3aT + т)и (п + 1, j + 1) + т и (п, j + 3)

- 2ти (п, j + 2) + (-4 h3 + г - 2 h3aт) и (п, j + 1) - 2 h2rF (п + 1, j + 2) + 2 h2r F (п + 1, j) - 2 h2r F (п, j + 2) + 2 h2rF (п, j)], [п, j], [uxx,ux,u,F]]

> collect(\%[1,1]/(4*tau*h**3),[tau,h,sigma]);

(-2 и (п + 1, j + 2) - 2 и (п, j + 2) V

+

22 -1 F (п + 1, j + 3) + 2 F (п, j + 1) + 2 F (п + 1, j + 1) - 2 F (п, j + 3)

1 / 1 1 Н 1

+ ^ (-4«(пз) + 2и(п +^ +1) - 2й(п + ^+ 3)

- 4и(п +1,з) + (п,з +1) + 4й(п,з + 4) - 1и(п,з + 3) + 4и (п + 1,, + 4)) + и (п + ^ + 2)-и (п^' + 2)

В результате получим следующую разностную схему для уравнения (14), аналогичную схеме Кранка-Николсона для уравнения теплопроводности

и"+1 - и] ^Й!1 -и^Жи3^ -и3?-!) + т 2 й +

«+21 - 2и"+11 + 2и]+! - <+2) + (и?+2 - 2и?+1 + 2и"_ 1 - и]_2)

+ - xj+2 ^"j + l "Г" ^"j-1 aj-2)~r\.d1j+2 2uj+1 + 2uj-1 uj-2

+ 4h

uj+1 + uj

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

^+1 = ^fc+1 - v3k + vl = (Vk+1 - Wfc)(vl+1 + Wfc+1 Wfc + vl) + vl « Wfc+1 ■ 3v2k - 2г»3-

Количество итераций для достижения точности 10-12 на следующем временном слое, как правило, не превышало 3. Шаг по времени t брался равным половине шага по переменной г]. Программа расчета была написана на языке Python с использованием пакета SciPy [12].

Результаты проведенного компьютерного моделирования представлены на рис. 2-4. Расчеты позволяют сделать следующие выводы. Наличие жидкости в оболочке приводит к существенному изменению характера распространения в ней продольных волн деформаций. Если в оболочке нет жидкости (эквивалентно условию а = 0), уединённая волна (кинк) движется, сохраняя свою первоначальную форму и скорость (см. рис. 2).

Наличие жидкости в оболочке из неорганических материалов ( а = 1) ведет к росту амплитуды волны (рис. 3). Таким образом, можно утверждать, что жидкость способствует постоянной дополнительной «подпитке» энергией (из источника первоначального возбуждения), обеспечивающей рост амплитуды.

0.20

0.15 0.10 0.05

0.00

S-

-0.05 -0.10 -0.15 -0.20

-°250 20 40 60 80 100

П

Рис. 2. График численного решения уравнения (14) с начальным условием

(15) при а = 0.0, к = 0.2 и для t = 0.0... 375.45

1.5

1.0

0.5

S- 0.0

-0.5

-1.0

150 20 40 60 80 100

П

Рис. 3. График численного решения уравнения (14) с начальным условием

(15) при а = 1.0, к = 0.2 и для г =0.0... 1.88

Наличие жидкости в оболочке из органических материалов (живые организмы, что соответствует а = -1) ведет к быстрому уменьшению амплитуды волны, т.е. к её затуханию (рис. 4). Для поддержки процесса распространения волны необходимо периодическое её возбуждение.

4. Заключение

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

— 4=0.00 ---4=75.09 - - 4=150.19 ...... 4=225.28 - 4=300.38 - 4=375.47 ■ ' / /

' ' / ' ' / ' / /

/ 7 1 / 1

/ 1 1

/ / 1 j 1 1

/ / I 1 1 1

I

— f=0.00 -- f =0.38 - -■- f =0.75 ...... f=1.13 - f=1.50 " - f=1.88

——-

- ' /

Рис. 4. График численного решения уравнения (14) с начальным условием

(15) при а = -1.0, к = 0.2 и для г =0.0... 1.88

Использование базиса Грёбнера для генерации разностной схемы при численном решении задачи Коши для нелинейного уравнения в частных производных третьего порядка по пространственной переменной позволило получить результат расчета без осцилляций, вызываемых численной реализацией. Численная схема также была протестирована на точном решении для а = 0 (см. рис. 2).

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

— рост амплитуды волны для неорганических материалов,

— падения амплитуды волны для живых организмов,

— отсутствие влияния жидкости для несжимаемых материалов, таких как резина.

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

Литература

1. Землянухин А. И., Могилевич Л. И. Нелинейные волны деформаций в цилиндрических оболочках // Изв. вузов. Прикладная нелинейная динамика. — 1995. — Т. 3, № 1. — С. 52-58. [Zemlyanukhin A. I., Mogilevich L. I. Nelineyjnihe volnih deformaciyj v cilindricheskikh obolochkakh // Izv. vuzov. Prikladnaya nelineyjnaya dinamika. — 1995. — T. 3, No 1. — S. 52-58. ]

2. Землянухин А. И., Могилевич Л. И. Нелинейные волны в цилиндрических оболочках: солитоны, симметрии, эволюция. — Саратов: Сарат. гос. техн. унт, 1999. — 132 с. [Zemlyanukhin A. I., Mogilevich L. I. Nelineyjnihe volnih v cilindricheskikh obolochkakh: solitonih, simmetrii, ehvolyuciya. — Saratov: Sarat. gos. tekhn. un-t, 1999. — 132 s. ]

3. Лойцянский Л. Г. Механика жидкости и газа. — М.: Дрофа, 2003. — 840 с. [Loyjcyanskiyj L. G. Mekhanika zhidkosti i gaza. — M.: Drofa, 2003. — 840 s. ]

4. Каузерер К. Нелинейная механика. — М.: Иностранная литература, 1961. — 240 с. [Kauzerer K. Nelineyjnaya mekhanika. — M.: Inostrannaya literatura, 1961. — 240 s. ]

5. Вольмир А. С. Оболочки в потоке жидкости и газа: задачи гидроупругости. — М.: Наука, 1979. — 320 с. [Voljmir A. S. Obolochki v potoke zhidkosti i gaza: zadachi gidrouprugosti. — M.: Nauka, 1979. — 320 s. ]

6. Блинков Ю. А., Мозжилкин В. В. Генерация разностных схем для уравнения Бюргерса построением базисов Грёбнера // Программирование. — 2006. — Т. 32, № 2. — С. 71-74. [Blinkov Yu. A., Mozzhilkin V. V. Generaciya raznostnihkh skhem dlya uravneniya Byurgersa postroeniem bazisov Gryobnera // Programmirovanie. — 2006. — T. 32, No 2. — S. 71-74. ]

7. Gerdt V. P., Blinkov Y. A., Mozzhilkin V. V. Grobner Bases and Generation of Difference Schemes for Partial Differential Equations // Symmetry, Integrability and Geometry: Methods and Applications. — 2006. — Т. 2. — С. 26. — http: //www.emis.de/journals/SIGMA/2006/Paper051/index.html.

8. Gerdt V. P., Blinkov Y. A. Involution and Difference Schemes for the Navier-Stokes Equations // Computer Algebra in Scientific Computing. — Springer Berlin / Heidelberg, 2009. — Т. 5743 из Lecture Notes in Computer Science. — С. 94-105.

9. LDA. — http://wwwb.math.rwth-aachen.de/Janet/lda.html.

10. Gerdt V. P., Robertz D. A Maple Package for Computing Grobner Bases for Linear Recurrence Relations // Nuclear Instruments and Methods in Physics Research. — 2006. — Т. A559. — С. 215-219. — ISSN 0168-9002. — ArXiv:cs.SC/0509070.

11. Maple. — http://www.maplesoft.com/.

12. SciPy. — http://www.scipy.org/.

UDC 531.383:532.516

Mathematical and Computer Modeling of Non-linear Deformation Waves in Shell with Viscous Liquid Inside

Yu. A. Blinkov*, S. V. Ivanov*, L. I. Mogilevich

* Mechanics and Mathematics Department Saratov State University named after N.G. Chernyshevsky Saratov, Russia, 410012 t Volga Branch Moscow State University of Railway Transport Saratov, Russia, 410790

The present investigation is devoted to the analysis of non-linear deformation waves propagation in physically non-linear elastic cylinder shell, containing viscous incompressible liquid. Wave processes in elastic cylinder shell without interaction with liquid were investigated earlier on the basis of soliton theory. The presence of liquid needs development of a new mathematical model and computer modeling of the processes, taking place in the system.

Key words and phrases: cylinder shell, oscillations, deformation waves, hydroelasticity, viscous incompressible liquid, solitary wave, Grobner basis.