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

УДК 533.6

И. И. Липатов1'2, В. К. Фам1

1 Московский физико-технический институт (национальный исследовательский университет) 2Центральный аэрогидродинамический институт имени профессора Н. Е. Жуковского

Моделирование панельного флаттера в рамках асимптотической теории течений вязкого газа

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

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

1'2 1

1

2

Modeling panel flutter in the framework of the asymptotic

theory of viscous gas flows

The processes of strong local viscous inviscid interaction in the flow around the flexible surface are studied. Linear and nonlinear processes of flow interaction in the laminar boundary layer with an external supersonic flow are investigated.

Key words: asymptotic theory of viscous gas flows, plate oscillation, fluid structure interaction.

1. Введение

Рассматривается обтекание плоской пластины сверхзвуковым потоком вязкого газа. Декартова система координат связана с пластиной, ось ОХ направлена вдоль поверхности пластины, ось ОУ по нормали к поверхности. Предполагается, что на расстоянии I от передней кромки расположен гибкий участок, имеющий длину д,.

Вводятся следующие обозначения для координат, отсчитываемых вдоль поверхности пластины и по нормали к ней, времени, компонентов вектора скорости, плотности, давления, полной энтальпии, коэффициента вязкости: 1х, 1у, И/и^, и^и, и^у, ржр, Р^и^р, и200Н/2, соответственно. Индекс те отно-

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

В зависимости от соотношения геометрических параметров и числа Ые возможны различные режимы течения. При воздействии возмущения давления с амплитудой Ар ^ 1 меняется толщина вытеснения пограничного слоя, при этом основной вклад в это изменение вносит пристеночная область с малыми скоростями [1]. При нелинейном воздействии на пристеночное течение изменения скорости можно оценить, как Аи ~ и ~ Ар1/2. Эта оценка верна, если влияние сил вязкости здесь в первом приближении несущественно.

© Липатов И. И., Фам В. К., 2019

(с) Федеральное государственное автономное образовательное учреждение высшего образования

«Московский физико-технический институт (национальный исследовательский университет)», 2019

При нелинейных изменениях скорости толщина пристеночной области сдвигового течения также меняется в главном члене, что следует из условия сохранения расхода. Тогда из оценки для продольной скорости у ~ ей, е = И,е-1/2 следует и оценка для изменения толщины вытеснения Ау ~ еАр1/2. При этом существенно, что основная часть пограничного слоя с конечными скоростями дает существенно меньший вклад в суммарное изменение толщины вытеснения А5 ~ еАр, поскольку изменения скорости здесь при малых амплитудах давления линейны.

Если возмущения, вносимые во внешний сверхзвуковой поток можно оценить по формуле Аккерета, и эти возмущения имеют тот же порядок, что и исходное возмущение давления, инициирующее процессы взаимодействия, тогда Ар ~ еАр1/2/Ах. Оценка для величины давления непосредственно определяет и длину области возмущенного течения Ах ~ е/Ар1/2. Отсюда следует, что для всех малых возмущений давления длина области возмущенного течения превосходит толщину пограничного слоя, что и предопределяет возможность использования формулы Аккерета для определения индуцированного возмущения давления.

Предположение о влиянии вязкости в области нелинейных изменений скорости приводит к известным оценкам теории свободного взаимодействия [1]. Ниже рассмотрен режим, для которого влияние вязкости в первом приближении несущественно. При малых возмущениях такой режим реализуется, если амплитуда давления удовлетворяет неравенству е1/2 << Ар << 1.

В этих условиях в поле возмущенного течения можно выделить 4 характерные области. Первая область содержит струйки тока невязкого сверхзвукового течения, характерный поперечный размер этой области определяется длиной этой области и наклоном характеристик, тогда для конечных чисел Маха у1 ~ е/Ар1/2.

Область 2 - это основная часть пограничного слоя. На дне этой области расположена область нелинейных возмущений продольной скорости (область 3), в которой влияние вязкости в первом приближении несущественно. Для учета влияния вязкости необходимо ввести в рассмотрение область 4, поперечный размер которой оценивается из условия баланса сил вязкости и инерции у4 ~ ез/2/Ар1/2. Следует отметить, что возможность существования предполагаемой структуры течения зависит от существования безотрывного течения в локальном пограничном слое (область 4).

2. Формулировка краевой задачи

Решение в области 3 можно записать в следующем виде, основываясь на полученных выше оценках:

^ = х3/р1/2а1/2рАр1/2, (1)

у = р-1/2а-1Ар1/2 уз, (2)

г = а-1р-1 Ар-Ч3, (3)

и(х, у, г, е, Ар) = гко-1/2Ар1/2из(хз,уз, 13) + ..., (4)

(х, у, г, £, Ар) = гко-1/2а-1/2рАр3/2Уз(хз, уз^з) + ..., (5)

р(х,у,г,е, Ар) = ) + Аррз(хз,Ьз) + ..., (6)

Р = Р-ш + ..., (7)

где параметр а определяется из решения для течения в невозмущенном пограничном слое а = ди/дуи,. Следует отметить, что этот параметр по порядку величины равен 0(Яе1/2). Подстановка выражений (1) - (7) в систему уравнений Навье-Стокса в результате предель-

ного перехода Re ^ те, Ар ^ 0, ApRe1/4 ^ те приводит к системе уравнений:

du3 + + + =0 (8)

dt3 dx3 dy3 dx3 ,

^ + ^ = 0, (9)

охг оуз

Р- = 0 (Ю)

дуз

с краевыми условиями

х ^ -те :v3 = 0,u3 = у3.

Решение можно искать в виде и3 = у3 + Аз(x3,t3). Уравнения (8) - (10) тогда можно привести к виду

— + А— + — = 0 (11)

dt дх дх '

где нижний индекс «3» опущен.

A

с обратным знаком. Во внешнем потоке это изменение индуцирует возмущение давления дА

А'Р = -TT-

дх

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

д дА ^ dw

дх дх'

дА дА др

+ + vw + / = 0. 13

dt дх дх

К этой системе следует добавить уравнение, описывающее деформацию гибкого участка [2], чтобы окончательно получить замкнутую систему уравнений:

^д4w „т92w d2w . .

Ds^w -NW + CW + **•t) = ° (")

w

Кроме того, необходимо уравнение для кинематической связи параметров в области 3 и пластины:

dw . dw

w, А

^д4 w ^д 2w d2w дА dw

дх4 дх2 dt 2 dx dx ' dA д®А d2 А d 2w dw ^dw

dt dx dx2 dx2 dt dx с граничными условиями

-1/2 <x < 1/2, w(x, t = 0) = gi(x); w(x = ±1, t) = 0,

dw(x = ±1, t) n dw(x,t = 0) , . -=-= 0;-dt-= 92(X),

/х дА(х —)• ±оо, £) А(х, i = 0) = да (х); ( ^-^ = 0,

где D = Eh3/ (12(1 — v2)) - цилиндрическая жесткость. N = Щ + Nx. Nq - натяжение. Nx = e— foi OWi dx - натяжение за счет прогиба пластины. E - модуль Юнга.

2 (X \ UX J

h

Для простоты будем решать систему уравнений (16) для безразмерных параметров D, N, С. dAw

— производная четвертого порядка имеет большую диссипативность и стремит-

ох4

d2w

ся стабилизовать колебание, а производная второго порядка -7—^7- имеет дестабилизующий

о х2

эффект. Нелинейный член в уравнении Бюргерса А—— характеризует процесс передачи

о х

о w

энергии от длинных волн к коротким волнам, нелинейный член А—— соответствует взаи-

ох

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

Рис. 1. Схема пластины

3. Решение линейной задачи

Линеаризуя систему уравнений (16), получим более простую систему уравнений:

+ дА + ды ^

дх4 0 дх2 Ы2 дх дх ' /171

дА д2А д2ы ды [ 4

+ + =0.

dt дх2 дх2 dt

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

А ~А0ei(^t-kx);W - ei("t-kx\

где Ао есть комплексное число.

Это приводит (17) к виду

i -cc2 + DkA + Nk2 + A0ik - ik = 0, [-A0icj -k2 - Ak2 - ico = 0.

Напишем А в виде А = a + bi. где a, b есть вещественные числа, получим выражение для прогиба пластины:

k2 (2Ь+ (1 -a2 - b2)i)

u =

b2 + (a + 1)2

Для того чтобы А, и> не возрастали со временем, мнимая часть частоты ш должна быть меньше либо равняться нулю 1т(ш) < 0, это соответствует условию а2 + Ъ2 > 1.

После нескольких математических преобразований получим следующую систему уравнений:

'(V - {а2 + Ь2 - + ык - ь

Ь2 + (а + 1)' 4Ь (1 - а2 - Ь2У_ а - 1

Ь2 + (а + 1)2 )2

к3

(19)

С помощью пакета МаЙаЬ [3] показано, что для всех волновых чисел к всегда существуют растущие моды. При И = 1, N = 1, к = 1 максимальный инкремент нарастания равен 0.4188 и частота -штах = 0.4188 + 1.6939г. С другой стороны, с помошью спектрального метода [4| можно привести систему линейных дифференциальных уравнений в частных производных к системе обыкновенных дифференциальных уравнений. Найдя все ее собственные числа, также получим растущую моду с наибольшим инкрементом нарастания, который равен 0.4188. Совпадение результатов подтверждает достоверность обоих методов.

При увеличении волнового числа максимальный инкремент нарастания уменьшается по линейному закону, а частота возрастает по квадратичному закону. При И = 1, N = 1, к = 5 получается -штах = 0.1 + 25.6г, а при И = 1, N = 1, к = 10 получается -^тах = 0.05 + 100г. Из этого следует, что неустойчивые моды длинных волн нарастают быстрее, чем неустойчивые моды коротких волн.

4. Решение нелинейной задачи

Дальше рассмотрена нелинейная задача, область пластины -1/2 < х < 1/2, область расчетов -20 < х < 20. Расчет проводится с помощью разностной схемы третьего и четвертого порядков точности и метода Рунге Кутта второго порядка точности.

Исследуется влияние дискретизации по времени и пространству на расчет, показано, что значения йх = 1е-2, йЬ = 1е-5 обеспечивают необходимую точность решений, поэтому для дальнейших расчетов используются последние дискретизации по времени и пространству.

Рис. 2. Дискретизация по пространству

Рис. 3. Дискретизация по времени

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

тов [5] выясняется, что существуют 4 главные моды, что эти моды содержат более 99% энергии колебаний пластины.

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

Рис. 4. Преобразование Фурье колебания пласти- Рис. 5. Нормированные амплитуды нормальных ны при х = 0.1 главных мод

Дальше показано влияние цилиндрической жесткости на колебание пластины. При И = 10 нелинейный член при производной второго порядка играет несущественную роль, поэтому колебание является регулярным. При И = 10~2 член при производной четвертого порядка мал, поэтому нелинейный член при производной второго порядка играет главную роль, так что колебание станет более сложным. Преобразование Фурье, фазовый портрет, вид главных мод при х = 0.1 доказывают это рассуждение.

Рис. 6. Главные нормальные моды Рис. 7. Преобразование Фурье и их главных мод

Рис. 8. Преобразование Фурье D = 0.01

Рис. 9. Преобразование Фурье D = 10

При И = 1е — 3 поведение пластины станет еще более сложным, с помощью метода Бубнова-Галеркина можно привести уравнение колебания пластины (уравнение фон Кармана) к системе обыкновенных дифференциальных уравнений. Разложим решение для прогиба пластины по собственным формам колебаний пластины Wj (х) в вакууме с неизвестными амплитудами Аj (I):

оо

ы(х, г) = YtWj (х^ (г). j=1

Рис. 10. Фазовый портрет И = 0.01

Рис. 11. Фазовый портрет И =10

Рис. 12. Главные моды И = 0.01

Рис. 13. Главные моды И =10

Подставим последнее выражение для прогиба пластины в уравнение фон Кармана в вакууме:

6Б

В4ы . лт лт х д2ы ды

— + ") ^ + С-ёё + ^ = 0,

ды

(ж

где С = 1, М

Умножим последнее уравнение на Wn (х) и проинтегрируем по х от —1/2 до 1/2. В силу ортонормированности собственных функций получим следующее уравнение для п-й амплитуды [6]:

где

дА

т2

п + и20пАп + 12 И ^ атк^пАтА^ = 0,

т,к^=1

1/2 íд4Wj „ д^ дWj

а

дх4

-М

дх2 дх

Wndx = ш^б1-,

1/2 -1/2

ШЩпс!х — 5™,

§п _ СИМвол Кронекера,

а]П -

1 [1

(дЩ дШЛ

\ дх дх )

йх — —■

1

г1/2 /д2щ

у/211-1/2\ дх2

йх.

Н/2

у/21 ¡-1/2

Считаем, что |Ага| << |А^,п > 1 и атп << а11,т,п > 1, тогда в уравнении для А\ можно отбросить все члены, которые содержат амплитуды с индексом выше 1. Тогда уравнение для А1 запишется в виде

д А 2

-д^ +<А1 +Ка2иА1 — 0, (21)

где К — 120.

Это есть известное уравнение Дюффинга, будем искать приближенное решение в виде асимптотического разложения типа Пуанкаре [7]:

и — и0 + еи1 + е и2...,

о о о

ш^ — ш2 + ес1 + е2С2... .

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

и0 + ш2и0 — 0, ио(0) — а, и0(0) — 0, (22)

и\ + ш2и1 + с1и0 + и0,и1(0) — 0, и 1(0) — 0. (23)

Решение уравения (22): и0 — а cos(шt)

ио

cos(3х) — 4cos(х)3 — 3cos(х), получим

и\ + ш2и1 + а(с1 + 3А2/4) cos(шt) + а3 )/4 — 0.

Отбросим вековой член, тогда С1 — —3 а2/4. Для первого приближения

22 ш221 — ш2 + е с1,

ш — ^ ш01 — 3е а2/4. Тогда получим решение уравнение Дюффинга в первом приближении:

и — а cos(ш ^ + а3/32cos(3шt),

где ш — у/ш21 — 3еа2/4, е — 12^^.

ш , 3 ш

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

Из графика видно, что для первой моды

ш11 — 6.2, ш12 — 18.8, ш12 & 3ш11.

Наблюдается внутрений резонанс двух первых мод:

(22 = 2.1, (11 = 6.2, (21 = 10.4, (Ш23 = 14.6, (12 = 18.8, ¿24 = 23, ш22 : (11 : (21 : (23 : (12 : (24 ~ 1:3:5:7:9:11.

Рис. 14. Внутрений резонанс между двумя первыми модами колебания

Рис. 15. Колебание пластины D = 1, А = 0.1Ч х = 0.1

Рис. 16. Колебание пластины D = 1, А =10, х = 0.1

Рис. 17. Бифуркационная диаграмма для локальных амплитуд колебания

Дальше исследуется влияние амплитуды возмущения толщины вытеснения пограничного слоя. Ясно, что увеличении амплитуды возмущения толщины пограничного слоя приводит к заметному нелинейному эффекту.

При А = 10 наблюдается бифуркация локальных амплитуд колебаний пластины (рис. 17) за счет сильной нелинейности.

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

Получены следующие выводы.

1. Существует критерий выбора главных мод колебания, заключающийся в том, что такие моды содержат более 99% энергии колебания. Зная эти моды, можно управлять колебанием пластины.

2. Наблюдается внутрений резонанс 1:3и1:3:5:7:9:11 между двумя первыми модами за счет нелинейности колебания пластины.

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

Литература

1. Нейланд В.Я., Боголепов В.В., Дудин Г.Н.; Липатов И.И. Асимптотическая теория сверхзвуковых течений вязкого газа. Москва : Физматлит, 2003.

2. Vedeneev V. V. Panel flutter at low supersonic speeds // Journal of Fluids and Structures. 2012. V. 29. P. 79-96.

3. https: / / www.mathworks.com

4. Trefethen L.N. Spectral Methods in MATLAB. SIAM, Philadelphia, 2000.

5. https://en.wikipedia.org/wiki/Principle component analysis

6. Веденеев В.В. Нелинейный высокочастотный флаттер пластины // Изв. РАН. МЖГ. 2007. № 5. С. 197-208.

7. Не J.H. Some Asymptotic Methods for Strongly Nonlinear Equations // International Journal of Modern Physics B. 2006. V. 20. P. 1141-1199.

References

1. Neiland V.Ya., Bogolepov V.V., Dudin G.N., Lipatov I.I. Asymptotic theory of supersonic viscous gas flows. Moscow : Fizmatlit, 2003. (in Russian).

2. Vedeneev V. V. Panel flutter at low supersonic speeds. Journal of Fluids and Structures. 2012. V. 29. P. 79-96.

3. https://www.mathworks.com

4. Trefethen L.N. Spectral Methods in MATLAB. SIAM, Philadelphia, 2000

5. https://en.wikipedia.org/wiki/Principle component analysis

6. Vedeneev V. V. Nonlinear high-frequency flutter of a plate. Fluid Dynamics. 2007. V. 5. P. 197-208. (in Russian).

7. He J.H. Some Asymptotic Methods for Strongly Nonlinear Equations. International Journal of Modern Physics B. 2006. V. 20. P. 1141-1199.

Поступим в редакцию 19.03.2019