ФИЗИКА, РАДИОТЕХНИКА И ЭЛЕКТРОНИКА
УДК 538.911
С.Г. Гестрин, Е.К. Сергеева ЧИСЛЕННЫЙ МЕТОД ИССЛЕДОВАНИЯ ДИФФЕРЕНЦИАЛЬНОГО УРАВНЕНИЯ РЭЛЕЯ
Предложен новый численный метод, основанный на совместном использовании метода Рунге-Кутта и метода бисекции, алгоритм которого основан на поэтапной корректировке одного из граничных условий, позволяющий найти решение уравнения Рэлея на отрезке, содержащем сингулярность.
Уравнение Рэлея, метод Рунге - Кутта, метод бисекции, ветровая неустойчивость, панельный флаттер
S.G. Gestrin, E.K. Sergeeva NUMERICAL METHOD FOR THE RESEARCH OF THE RAYLEIGH DIFFERENTIAL EQUATION
We propose a new numerical method based on the joint use of the Runge - Kutta and bisection method. This algorithm is based on the incremental adjustment of one of the boundary conditions to find the solution for the Rayleigh equation on the interval containing singularity.
Rayleigh equation, Runge-Kutta method, bisection method, wind instability, panel flutter
Одна из основных задач теории гидродинамической неустойчивости - это задача о генерации волновых возмущений сдвиговыми потоками [1]. Наиболее сложным в этой связи является вопрос о взаимодействии волн различной природы с течением, в котором скорость меняется непрерывным образом.
В [2] было показано, что нарастание волн на поверхности глубокой воды связано с их резонансным взаимодействием с воздушным течением над поверхностью. Резонансное усиление поверхностных волн происходит в критическом слое сдвигового потока, где его скорость U 0 (yc) близка к фазовой скорости поверхностной волны. Данный механизм называют ветровой неустойчивостью.
В [3, 4] было показано, что механизм ветровой неустойчивости может приводить к усилению колебаний пластины при ее взаимодействии со сдвиговым газодинамическим течением. Данный эффект называют панельным флаттером.
Рассмотрим несжимаемую идеальную жидкость, находящуюся между колеблющейся пластиной, поверхности которой соответствует y = y1 = 0 , и твердой стенкой y = y2 . Жидкость движется
со скоростью V0 = (U0 (y),0). Амплитуда колебаний функции тока у при этом, как известно, удовлетворяет уравнению Рэлея [1 - 4]:
d 2у С U '(y) ^
+ k
У = 0. (1)
dy2 ^ U ( у )-ю/ k
Заметим, что уравнение (1) играет центральную роль в теории неустойчивости сдвиговых течений. Оно содержит особенность в точке yC, где скорость потока совпадает с фазовой скоростью волновых возмущений: U0 (yC) = $/ к .
Для построения решения уравнения Рэлея применим численный метод, в основе которого лежат метод Рунге-Кутта и метод бисекции [5]. Запишем уравнение (1) в приближенной форме, разлагая U 0 (у ) в ряд Тейлора вблизи точки Уc :
d У f Uо(yC ) , + k 2
л
у = 0. (2)
Лу2 IU0(Уc)(у - Уc)
Представляя уравнение (2) в виде эквивалентной системы двух уравнений первого порядка для функций w1 = у/(у)/У(0), w2 = (у(у)/^(0)) , получим
w/1 = w2,
4 в
w 2 = wl
f в 2 л U0(Ус ) f wl(0) =1^
і У — y,
U 0 (yc )’ і w2(0)
wl
(З)
Граничное условие wl(0) = 1 является очевидным, в то время как w2(0) не определено. В дальнейшем будем искать решение системы (З) на отрезке [yl, y2], где у2=0. Потребуем также выполнения условия wl(yl) = 0. Физически оно означает, что на расстоянии у1 от колеблющейся пластины находится твердая стенка, на поверхности которой нормальная к ней компонента скорости vy (yl)=—dy/dx(yl) = 0.
Значение w2(0) должно быть подобрано таким образом, чтобы wl(yl) = 0. Для примера приведем результаты в случае yl = —Т , k = 1,5, в = —0,05, yc = —1. Для начала в качестве пробного значения возьмем w2 (0) = l. При отсутствии резонансного слагаемого в уравнении Рэлея в — 0, из условия w2(0) = 1 следует, что решение представляет собой при у < 0 одну лишь затухающую с удалением от пластины экспоненту у(у) = y(0)exp(ky). С появлением резонансного слагаемого в ^ 0 появляется примесь также нарастающего решения.
На рис. l приведен результат поэтапной корректировки граничного условия w2(0) . В виде сплошной линии представлена зависимость wl(у), линия из точек - w 2(y) , пунктирная линия -
w(у) = exp(l,5y), w^ — у . Для расчетов использовался метод Рунге-Кутта четвертого порядка с постоянным шагом на сетке из З0000 равноотстоящих узлов. В качестве первого пробного значения полагаем w2(0) =l и получаем решение системы (З) по методу Рунге-Кутта (см. рис. l а).
Обращаем внимание на знак решения в точке yl = —Т . Для изменения w2(0) задаем некоторый шаг h, например, Aw(0) = h = 0,5 . Снова ищем решение (З) уже для нового значения w2(0) = 1 + 0,5 = 1,5 (см. рис. 1 б). Если знак решения в точке yl = —Т не изменился, то сдвигаемся еще на один шаг h вправо, при этом новое значение w2(0) = 2, если же знак изменился, то делим шаг h/ 2 пополам, сдвигаемся влево и ищем решение для w2(0) = 1,5 — 0,25 = 1,25 (см. рис. 1 в). При этом новый шаг Aw(0) = 0,25. Продолжаем эту процедуру, смещаясь при изменении знака решения, в точке yl = —Т на полшага назад. В случае если знак решения не изменяется, то на полшага вперед (см. рис. 1 г, д) и т.д. На рис. 1 е приведен результат на 16-м этапе.
Видно, что в точке yl = —Т значение функции wl(— Т) = —0,0625 уже совсем мало отличается от 0. Продолжая выполнение данного алгоритма, можно добиться того, чтобы условие wl(— Т) = 0 выполнялось с любой наперед заданной степенью точности |wl( (yl) < Є .
На рис. 2 представлены окончательные результаты расчетов для в = —0,З;—0,Т . Видно, что с увеличением значимости резонансного слагаемого (ростом |в|) отличие wj (у) от exp(ky) становится более заметным, чем на рис. 1 е. Начиная с определенных значений в, поведение решения меняется качественным образом. Если при малых в решение всюду убывает с удалением от у 2 = 0, то
при больших в наблюдается рост вплоть до резонансного слоя, а лишь затем начинается спад (см. рис. 2 б).
78
c
Рис. 1. Расчетные кривые при в = -0,05: а - w2(0) = 1; б - w2(0) = 1,5; в - w2(0) = 1,25 ; г - w2(0) = 1,375 ; д - нф) = 1,4375 ; е - w2(0) = 1,4700938
б
а
в
г
е
д
-4 -2
а б
Рис. 2. Расчетные кривые при к = 1: а - р = -0,3, ш2(0) = 0,6437 ; б) р = -0,7 ш2(0) = -1,0368
Таким образом, создан численный метод, основанный на методе Рунге - Кутта и методе бисекции, позволяющий найти решение уравнение Рэлея на отрезке, содержащем сингулярность. Показано, что поведение решения существенно зависит от кривизны профиля скорости в критическом слое.
ЛИТЕРАТУРА
1. Степанянц Ю.А. Распространение волн в сдвиговых гидродинамических течениях / Ю.А. Степанянц, А.Л. Фабрикант // Успехи физ. наук. 1989. Т. 159. Вып. 1. С. 83-123.
2. Miles J.W. On the generation of surface waves by shear flows // J.W. Miles. Fluid Mech. 1957. V. 3. P. 185-204.
3. Гестрин С.Г. Ветровая неустойчивость и резонансное взаимодействие упругих колебаний тонкой пластинки со сверхзвуковым газодинамическим потоком / С.Г. Гестрин, Е.К. Сергеева // Известия вузов. Физика. 2011. № 3. С. 89-94.
4. Гестрин С.Г. Резонансное взаимодействие упругих колебаний тонкого стержня со сдвиговым течением «мелкой воды» / С.Г. Гестрин, А.Н. Сальников, Е.К. Сергеева // Известия вузов. Физика. 2010. № 1. С. 28-33.
5. Плис А.И. MATHCAD математический практикум / А.И. Плис, Н.А. Сливина // М.: Финансы и статистика, 1999. 656 с.
Гестрин Сергей Геннадьевич -
доктор физико-математических наук, профессор кафедры «Физика»
Саратовского государственного технического университета имени Гагарина Ю.А.
Sergey G. Gestrin -
Dr. Sc., Professor
Department of Physics
Gagarin Saratov State Technical University
Сергеева Елена Константиновна -
аспирант кафедры «Физика»
Саратовского государственного технического университета имени Гагарина Ю.А.
Elena K. Sergeeva -
Postgraduate
Department of Physics
Gagarin Saratov State Technical University
Статья поступила в редакцию 24.10.11, принята к опубликованию 01.12.11