CC BY
55
МАТЕМАТИКА И МЕХАНИКА
УДК 681.03.06:531.383:532.516
А.Ю. Блинкова
МОДЕЛИРОВАНИЕ НЕЛИНЕЙНЫХ ВОЛН ДЕФОРМАЦИЙ В ФИЗИЧЕСКИ ЛИНЕЙНЫХ ВЯЗКОУПРУГИХ ЦИЛИНДРИЧЕСКИХ ОБОЛОЧКАХ,
СОДЕРЖАЩИХ ВЯЗКУЮ НЕСЖИМАЕМУЮ ЖИДКОСТЬ
Настоящее исследование посвящено анализу распространения нелинейных продольных волн деформаций в цилиндрической оболочке, содержащих вязкую несжимаемую жидкость внутри оболочки. Физические свойства оболочки определяются уравнениями линейной теории вязкоупругости. Наличие жидкости потребовало разработки новой математической модели и компьютерного моделирования процессов, происходящих в рассматриваемой системе.
Цилиндрическая оболочка, колебания, волны деформации, гидроупругость, вязкая несжимаемая жидкость, солитон
A.Yu. Blinkova
NONLINEAR WAVES DEFORMATION MODELLING IN PHYSICALLY LINEAR VISCOELASTIC CYLINDRICAL SHELL WITH VISCOUS INCOMPRESSIBLE FLUID INSIDE
The present study is devoted to the analysis of nonlinear deformation of longitudinal waves in a cylindrical shell containing viscous incompressible fluid inside.
The physical properties of the shell are determined by the equations of the linear theory of viscoelasticity. The presence of fluid required working out of a new mathematical model and computer simulation of the processes occurring in the system.
Cylinder shell, oscillations, deformation waves, hydroelasticity, viscous incompressible liquid, solitary wave
1. Волновые процессы в вязкоупругих и нелинейно вязкоупругих оболочках, не взаимодействующих с вязкой жидкостью, рассмотрены в [1-3].
Получим, уравнения динамики с учётом наличия вязкой несжимаемой жидкости в
цилиндрической оболочке с помощью асимптотических методов для решения связанной задачи гидроупругости с соответствующими граничными условиями.
Рассмотрим бесконечно длинную вязкоупругую цилиндрическую оболочку, внутри которой находится вязкая несжимаемая жидкость.
Уравнения движения вязкой несжимаемой жидкости и уравнение неразрывности в цилиндрической системе координат r, x записываются в случае осесимметричного течения в виде [4]
+ grad1 V2 + rotV xV +—grad ■ p = -v rot rotV, m
dt 2 p (1)
divV = 0.
На границе с оболочкой выполняются условия прилипания жидкости
v v dU и ш пл
Vr = —— ,Vx = — при r = R1 - W, (2)
dt dt
Здесь t - время; Vr,Vx - проекции вектора скорости жидкости на оси цилиндрической системы координат; p - давление; p - плотность; V - кинематический коэффициент вязкости; U -продольное упругое перемещение оболочек по оси x ; W - прогиб, положительный к центру кривизны оболочки; R1 - внутренний радиус оболочки.
В случае осевой симметрии, используя гипотезу Кирхгофа-Лява, имеем связь между компонентами деформаций £х, £у перемещениями [5]
ди д2Ш 1 .ди дШ 2 1 дW^2 1 Т17 ^
£ = аТ-чхг+2("дХ-+ ^ £ , (3)
где Я - радиус срединой поверхности оболочки, г - расстояние от нее. Связь между компонентами напряжений Ох, ау и деформаций зададим уравнениями линейной теории вязкоупругости [6],
учитывающей линейную упругость объёмных деформаций
Е , Е
а,
£ +№у) —,-----а\‘ е-р{(-Т)вхйт,
у 1 _|_ //
х ~ , 2 4 х ' г^О у> ,
1 -Мо 1 + Мо (4)
^ =тЫ(еу +ДА)
Здесь Е - модуль Юнга, /И0 - коэффициент Пуассона материала оболочек (считая их одинаковыми), ( - время; а, в - параметры вязкоупругости; ех, еу - компоненты девиатора деформаций
2 4 / 2 2 21 2 1
£и =3(£х + £у -£х£у); ех =3 £х - 3£у , еу =3 £у - 3 £х. (5)
Разлагая функции ех, еу в ряд Тейлора по степеням (( - Т), при условии в( >> 1 сохраняем два члена разложения из формул (4) получим приближенные уравнения состояния [1-3]
Е , г2 1
ах = “--------2(£х +Мо£у ) + РЬ £х -ъ£у ]
1 - Ао 3 3
Е , ч г2 1 п
ау = + А£х) + А^у -3£х]
(6)
где введен оператор р, такой, что
г Е / а д а \ { С74
Р/ =-л-----(^2^ -~Б)/ (7)
1 + А0 в д( в
Вычислим с использованием (6) усилия и моменты по формулам [5]
Д Д Д Д
= 110аА’Ny = I 10ау^,мх = 110ахгёг,Му = |\ауггйг (8)
h а М — I 2 а гЙг ЛМ — I 2
2 2 2
и подставим (8) в систему уравнений динамики оболочек [5]
д^х , д 2и д 2МХ 1 д дШ ... . д 2W
17 ' рА Шг = ~“х ’ ~дхг+~Яыу+дх (~ъ1ых) - =~ч- (9)
здесь h0 - толщина оболочки; чх, чп - напряжения, действующие со стороны жидкости на поверхность оболочки, снесенные на невозмущенную поверхность оболочки (Ш << Я)
Чх = [РК^ + д^)]г=я,дп = [-р + 2рУд^]г=к (10)
дг дх дг
2. Принимая за характерную длину волны деформации I, перейдем к безразмерным
переменным для исследования уравнений динамики оболочек (3)-(9)
— х I Е
Ш и = ити1. (* = -Т(, х* = -, с0=л—Т1-21. (11)
I I V Ро (1 -АО )
здесь с0 - скорость звука в материале оболочки.
Положим
= £ = 0(1), ^ = 0(£), а = 0(1),
1 Я в (12)
т = °(£), = °(£ш)- hТ = 0(£),
где £ << 1 - малый параметр задачи.
Применим метод двухмасштабных разложений, вводя независимые переменные в виде
2
с* & & &
д = х- с( , т = £(, (13)
где - - безразмерная неизвестная скорость волны, а зависимые переменные представлены в виде
разложения по малому параметру £:
и = ию + £ии +..., и3 = и30 + £<31 +. .. (14)
Подставляя (11), (13), (14) в уравнения (3-9) с учетом оценок (12), получим в нулевом
приближении по £ линейную систему уравнений, из которой следует связь
Я
иш ди10
I дд ’а 1
1/1 ча
ао + 3(1 -Мо) в
(1 -Мо)
а
в
(15)
и определяется безразмерная скорость волны
с2 = [1 -1(1 -а) |](1 -а2)- (16)
Из следующего приближения по £, учитывая (15) и (16), находятся уравнения для определения и10:
2 - Я2 д Чо + -иж ди1Т д Чо + а1 ^ _ ,.„ +
ддэт 21 £ дд 21£ дд дд
1 ^ ч а с0 ^ ^ д3и
3
-(1 -Ао)Т°(1 + А1 +М\) ез , 20
в £ дд ЕитРХ-ъ2-
ПО
-—^[Чх -А Я%]
(17)
е дд
3. Для определения правой части уравнения (17) введем безразмерные переменные и параметры
-0 ,, со * г * С0 + * х Р^с01^т Г}
V = w — V , V = w — V , г = —, ( = — (, х = —, р = ------------“ т Р;
г т 1 г ’ х т^ х' т~» ’ 1 ’ ^ ’
/
Я
/
Я 1
Щ = — = 0(£2), Л = —,
/
Я13
I
Я
(18)
/ << 1, Л << 1.
Подставляя (18) в уравнения гидродинамики (1) и граничные условия (2), представим безразмерные скорость и давление в виде разложения по малому параметру Л :
V 0 +Лу„ +Р = Р0 +ЛР1 +...
(19)
В нулевом приближении по Щ (Щ ~ 0 - гидравлическая теория смазки), считая
//\Я1с0/у\<< 1 (ползущие течения [6, 7]), и в нулевом приближении по Л получаем уравнения гидродинамики (классические уравнения гидродинамической теории смазки)
дР0 _ 0 дР0
дг
= 1 д * дvХ0.
* —. * V —ч * /? *
дх г дг дг г дг
1 д * 0 дvХ
(г V;0)+ ^ = 0
дх
и граничные условия
дУ
Э^0
г тг~* = 0, г ^ дг дг
0 ди
=0
*30
Я, ди
, V
д(* х
10
^1 д(
Из решения задачи (20), (21) следует, что
1 итЯ1 ди10
при г* = 0,
* 7
при г =1.
(20)
(21)
Ри
дг
16|
•ди
30
2 д( *
д( *
йх
йх\
дг*
= 1 * дР0
2 дх*
(22)
С принятой точностью по £,Щ,Л из (5) найдем
и
V
V
х
г
0
*
wmc0 дvx . wтc0 I п0
чх = ру т2° х* I * 1, чп = -ру-^°—р ,
Я дг г =1 Чп ^2 Я1
и выражение в квадратных скобках правой части (22) имеет вид
Я дапл wmc0г/дv“ Я I дР\ wmc0 1 дР° „ Я п
4 "А 7э£]=^-Яг1(аг*)r•_l +А тяг а? = 2[1+2А я;] (23)
Учитывая, что были введены переменные (13), (14) и имея соотношения (15) (16), имеем
ро = 8сиЛ[2А1Я- - 1]и10, с =1[1 - 2(1 -Ао)^](1 -а2) . (24)
Я1 V 3 в
Следовательно, в правой части уравнения (17) остается выражение
2 --^-[1 - (2А1Я )2] (25)
АЛ Я1с0£ Я1 дд
с принятой точностью по Щ, £ положим Я1 ~ Я.
Подставляя (23) в уравнение (17), окончательно получим
дЧо + дЧо + 1(Я)2 2 с дЧо -
дддт 1£ 2 дд дд2 £ г 1 2 дд4
- 3 а Г(1 “Ао)(1+А+А12) % - 2[1 - (2д)2] ^ = 0.
3 в 1£ дд -,дя1с0£ дд
При отсутствии жидкости (— = 0) последнее слагаемое выпадает и уравнение превращается в
ди
уравнения Кортевега-де Вриза-Бюргерса, (КдВ-Б) для ---, имеющие точное частное решение. В
дд
зависимости от физических параметров величина А может быть больше 1, меньше 1 или равна 1.
2 2 2
Последний случай эквивалентен отсутствию жидкости, но означает, что она не влияет на волну деформации.
Легко видеть, что замена
ди1и _
— -3р, ц = c1д, ( = с2
(26)
дд
позволяет записать уравнение (26) в виде
др , др д р д 2р
— + 6р—+—Чг-а2 —г
д( дц дп дц
10 = с3р, ц = с^, ( = с2Т (27)
——+ 6р——+ ——3—аг-—2—ар = 0. (28)
1 1 1 п 1
здесь а = 1 при А1 < —, а = -1 при А1 > 2 и а = 0 при А1 = 2.
Постоянные с3, с15 с2 определяются формулами
с2=2а[1 - (2а)2] —~ , с1 =
ро!^о Я1со£
I \2 2
с2Л я) А
1/3
, с2 21£
, с3 = 6—-------------------
с1 сит
а2 = с1-—^77(1 -Ао)(1+А1 + А). (29)
При этом вводятся обозначения
1 а с0 3 в2 ~£
4. Запишем уравнение (28) в интегральной форме
£(-3р2 - Рцц + а2Рц )й( + рйц - Царёхйц = 0 (30)
для любой области ^. Для перехода к дискретной формулировке сопоставим и* = р((п ,Ц ■) и выберем в качестве базового контур, показанный на рис. 2.
Рис. 1. Базовый контур для уравнения (28)
Добавим интегральные соотношения
г
Jn
]uvdn - u(tn) — u(t,n}),
j
J^j+1unn dn - un(t’ Піі) — un(t’ Vj) •
(31)
Используя для интегрирования по времени и по четным производным по П формулу трапеций, а по нечетным производным по Т] формулу среднего значения, и полагая їп+1 — їп = Т, П+1 — П^ = Ь, перепишем соотношения (30), (31) в виде
(—З(
а, 2» 2
3\u j і u j
n 1І
u j 12 — u j 12
n 12)—(
n n1l n n1l I
\unn j 1 unn j — unn j 12 — unn j 12>
(n , n1l n
u 1 u — un
■ n] n] n] 12
— u4j 1l2))--1 (u” 11 — ■ 2h — 1 u» 11) ■hT - О
h
(32)
nj 11
2
nn j 11— uJ,
n r\ l n n
unn . , • 2h = un. - un ..
nn]+1 n]+2 V]
В результате разностная схема для уравнения (28) определяется как условие совместности для данной системы разностных уравнений (32). Таким образом, получается разностная схема [9, 10, 11], автоматически обеспечивающая выполнение интегральных законов сохранения по областям, составленным из базовых конечных объемов.
Для построения разностной схемой воспользуемся приведенной ниже программой, написанной на языке системы компьютерной алгебры Singular [12]. ring r = (0, h, tau, sigma2, sigma), (Tx, Tt), (c, dp);
// u_xx, u_x, u, 3uA2
vector eq1 = [-(1+Tt-TxA2-Tt*TxA2)*tau/2, sigma2*(1+Tt-TxA2-Tt*TxA2)*tau/2,
(1-Tt)*2*h - sigma*(1+Tt)*h*tau,-(1+Tt-TxA2-Tt*TxA2)*tau/2]; vector eq2 = [0, (Tx+1)*h/2,-(Tx-1), 0]; vector eq3 = [Tx*2*h,-(TxA2-1), 0, 0]; module m = eq1, eq2, eq3; std(m);
В первой строке программы описан полиномиальный модуль с переменными Tx, Tt с исключающим по позиции упорядочением над кольцом рациональных чисел с параметрами h, tau, sigma 2, sigma. Как видно из следующего комментария программы, первой позиции
соответствует функция uxx, а затем по порядку ux, u,3u2. Переменные Tx, Tt соответствуют операторам сдвига по переменным TJ и t. За счет выбора исключающим по позиции упорядочения
нелинейная часть 3u 2 не будет входить в лидирующие мономы системы при построении базиса Грёбнера командой std(m). В приведенном ниже результате вычислений первый элемент базиса Грёбнера представляет собой искомую разностную схему для уравнения (28), аналогичную схеме Кранка-Николсона для уравнения теплопроводности.
_[1]=[0,0,(tau)*TxA4*Tt+(tau)*TxA4+(-2*h*tau*sigma2-2*tau)*TxA3*Tt+ (-2*h*tau*sigma2-2*tau)*TxA3 +(4*h*tau*sigma2)*TxA2*Tt+ (4*h*tau*sigma2)*TxA2+(-2*hA3*tau*sigma-4*hA3-2*h*tau*sigma2+
2*tau)*Tx*Tt +(-2*hA3*tau*sigma+4*hA3-2*h*tau*sigma2+2*tau)*Tx+
n
(-tau)*Tt+(-tau), (hл2*tau)*TxЛ3*Tt+(hл2*tau)*TxЛ3+(-hл2*tau)*Tx*Tt+ (-^2*tau)*Tx]
_[2]=[0,(h)*Tx+(h),-2*Tx+2]
_[3]=[(-2*h*tau)*Tt+(-2*h*tau),(tau)*TxЛ3*Tt+(tau)*TxЛ3
+(-2*h*tau*sigma2)*TxЛ2*Tt+ (-2*h*tau*sigma2)*TxЛ2+(-tau)*Tx*Tt+ (-tau)*Tx+(2*h*tau*sigma2)*Tt+(2*h*tau*sigma2), (-4*^2*tau*sigma-8*^2)*Tt+(-4*^2*tau*sigma+8*^2), (2*h*tau)*TxЛ2*Tt+(2*h*tau)*TxЛ2+(-2*h*tau)*Tt+(-2*h*tau)] _[4]=[(2*h)*Tx,-TxЛ2+1]
Перепишем полученную разностную схему в обычных обозначениях
un+1 un I 2n1i 2п11ч . I 2n 2n \
u ■ - u ■ + 3 (u j+i - и j-l) + (u j+1 - и j-l) +
т 4Н
+ (и% - 2И"‘ + 2И" - И")+(И;3 - 2И;,+2И;_, - «■_,> _
+ 4^ (33)
/ П+1 /-> П+1 . П+1 \ . / П О П . П \ п+1 . п
о и+1 - 2и; + и-1) + (и;+1- 2М + и1-1) о и} + и} 0
- о9----------------------------------------------------------«-о-— и.
2 2Н1 2
Для определения свойств полученной разностной схемы (33) построим ее 1-е
дифференциальное приближение [13]
и( + 6иип + и^пп - о^Мщ - ои — (1/12)т2о3и + (- 108т2ииЩ + 36т2иоищ - (7/2)т2ио2)Мщ +
+ ((-99т2Мщ + 6т2(о + 602и))ищщ + (54т2о1иП1 + (-162т2и2 - 3Н2 - (17/2)т2о2о)Мщ +
+ (1/4)т2о(о2о + 72м2 )))ищщ + (27/2т2о:иппп + (-т2 (99и + 8о2 )Мщщ - 81т2иЩ +
+ 3/2т2(36о2и + 7о)ищ - Н2и - 1/4т2о2 - 9/2т2о2ио- 18т2(и)3))ипщ +
+ (-21т Мщщп + 21т о2Мщп - (9/2)т (о2 + 14и)Мщ+ 9т о2и +
+ (1/4)т2о22о+ (9/2)т2ио+ (1/12)^2о2)ипппЩ + (-(3/2)т2о22и - (1/4)й2 -
- (1/2)т2о2о - 9т2м2 + (21/2)т2о2Мщ - (27/2^^)^^ +
+ (-6т2ип + (1/12)т2(36о2и + 3о + о\))ищщЩ -
- (1/4)т (6и + о2 )ищпппппп + 1/4т о2МппПППППП - (1/12)т2иппппппппп Из построенного 1-го дифференциального приближения видно, что при достаточно малом шаге по времени т разностная схема (33) хорошо аппроксимирует уравнение (28). 5. Полученная неявная разностная схема имеет квадратичную нелинейность для следующего временного слоя. При построении ее решения методом простой итерации использована следующая линеаризация:
ук+1 — у2+1 - У1 + у1 — у +1 - ук )у+1 + у) + V2 ® Ук+1 • Ък - У2к.
В качестве начального условия при решении задачи Коши для уравнения (28) можно выбрать следующее точное решение уравнения при о — 0.
5 О 1 2 1 2 1
(р =-----+ — a22-------a^tanh(9)-apanh (9), 9 = — J2х -°о (34)
3 j2 50 25 50 10
Шаг по времени t брался равным половине шага по переменной Т). Программа расчета написана на языке Python с использованием пакета SciPy [14].
Результаты проведенного компьютерного моделирования представлены на рис. 2-4. Наличие жидкости в оболочке приводит к существенному изменению характера распространения в ней продольных волн деформаций. Если в оболочке нет жидкости (эквивалентно условию J = 0), уединенная волна (имеет структуру ударной волны) движется, сохраняя свою первоначальную форму и скорость (см. рис. 2).
Наличие жидкости в оболочке при j = 1 ведет к росту амплитуды волны (см. рис. 3). Таким
образом, можно утверждать, что при при и < — жидкость способствует постоянной дополнительной
2
«подпитке» энергией (из источника первоначального возбуждения), обеспечивающей рост
амплитуды.
Рис. 2. График численного решения уравнения (28) с начальным условием (7) при о — 0.0, о2 —1.0, о — 0.5 и для г — 0.0... 7.52
Наличие жидкости в оболочке при о — -1 ведет к быстрому уменьшению амплитуды волны,
то есть к её затуханию (см. рис. 4). Для поддержки процесса распространения волны при и > —
2
необходимо периодическое её возбуждение.
1------------------------------!-----------------------------Г
Рис. 3. График численного решения уравнения (28) с начальным условием (7) при о — 1.0, о2 —1.0, о— 0.5 и для г — 0.0.1.00
Заключение
Проведенное моделирование с использованием компьютерной алгебры позволило выявить особенности поведения волн в физически линейных вязкоупругих цилиндрических оболочках, содержащей вязкую несжимаемую жидкость.
Использование базиса Грёбнера для генерации разностной схемы при численном решении задачи Коши для нелинейного уравнения в частных производных третьего порядка по пространственной переменной позволило получить результат расчета без осцилляций, вызываемых численной реализацией. Численная схема также была протестирована на точном решении для о — 0 (см. рис. 2).
Рис. 4. График численного решения уравнения (28) с начальным условием (7) при a = -1.0, j2 = 1.0, (О= 0.5 и для t = 0.0...1.00
Полученный расчет показал влияние вязкой несжимаемой жидкости на поведение нелинейной волны деформации в оболочке в зависимости от величины и , характеризующей материал оболочки:
1 1 рост амплитуды волны для и < , падения амплитуды волны для и1 > ^, отсутствие влияния
жидкости для и = —. За счет рассеяния энергии в вязкоупругом материале оболочки происходит
сглаживание профиля волны деформации (см. рис. 3, 4).
В заключение приношу благодарность профессору Л.И. Могилевичу за постановку задачи и внимание к работе.
Работа выполнена при финансовой поддержке гранта РФФИ № 13-01-00049а и Гранта Президента РФ № МД-1025.2012.8
ЛИТЕРАТУРА
1. Землянухин А.И. Нелинейные волны в цилиндрических оболочках: солитоны, симметрии, эволюция / А.И. Землянухин, Л.И. Могилевич. Саратов: СГТУ, 1999. 132 с.
2. Аршинов Г.А. Двумерные уединенные волны в нелинейной вязкоупругой деформируемой среде / Г.А. Аршинов, А.И. Землянухин, Л.И. Могилевич // РАН. Акустический журнал. 2000. Т. 46. № 1. С. 116-117
3. Аршинов Г.А. Статические и динамические задачи вязкоупругости / Г.А. Аршинов, Л.И. Могилевич. Саратов: Сарат. гос. агр. ун-т, 2000. 152 с
4. Лойцянский Л.Г. Механика жидкости и газа / Л.Г. Лойцянский. М.: Дрофа, 2003. 840 с.
5. Вольмир А.С. Нелинейная динамика пластинок и оболочек / А.С. Вольмир. М.: Наука, 1972. 432 с.
6. Москвитин В.В. Сопротивление вызко-упругих материалов / В.В. Москвитин. М.: Наука, 1972. 328 с.
7. Чивилихин С.А. Динамика скручивающихся нанотрубок в вязкой жидкости / С.А. Чивилихин, И.Ю. Попов, В.В. Гусаров // Доклады РАН. 2007. Т. 412. № 2. С. 201-203.
8. Солитоны в стенке нанотрубки и стоксово течение в ней / Ю.И. Попов, О.А. Розыгина, С.А. Чивилихин, В.В. Гусаров // Письма в ЖТФ. 2010. Т. 36. Вып. 18. С. 42-54.
9. Блинков Ю.А. Генерация разностных схем для уравнения Бюргерса построением базисов Грёбнера / Ю.А. Блинков, В.В. Мозжилкин // Программирование. 2006. Т. 32. № 2. С. 71-74.
10. Gerdt V.P. Grobner bases and generation of difference schemes for partial differential equations / V.P. Gerdt, Yu.A. Blinkov, V.V. Mozzhilkin // Symmetry, Integrability and Geometry: Methods and Applications. 2006. Vol. 2. P. 26. http://www.emis.de/journals/SIGMA/2006/Paper051/index.html.
11. Gerdt V.P. Involution and difference schemes for the Navier-Stokes equations / V.P. Gerdt, Yu.A. Blinkov // Computer Algebra in Scientific Computing. Springer. Berlin: Heidelberg, 2009. Vol. 5743 of Lecture Notes in Computer Science. P. 94-105.
12. Singular. http://www.singular.uni-kl.de/
13. Шокин Ю.И. Метод дифференциального приближения. Применение к газовой динамике / Ю.И. Шокин, Н.Н. Яненко. Новосибирск: Наука, 1985. 363 c.
14. SciPy. http://www.scipy.org/
Блинкова Анастасия Юрьевна -
аспирант кафедры «Теплогазоснабжение, вентиляция, водообеспечение и прикладная гидрогазодинамика» Саратовского государственного технического университета имени Гагарина Ю.А.
Статья поступила в редакцию 03.08.12, принята к опубликованию 06.11.12
Anastasiya Yu. Blinkova -
Postgraduate
Gas Supply, Ventilation, Water Supply
and Applied Fluid Dynamics
Gagarin Saratov State Technical University
