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

УДК 539.3:533.5:517.9

П. А. ВЕЛЬМИСОВ, С. В. КИРЕЕВ

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

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

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

г а о V2 V РР Г>м/(4) +N»' + <1*/+ /(}*)-№ |(У)2^х = 0> а = ,оРо , М =—, 0 =-=—7-, П=Е1 (1)

О л/л^-1 я 2£(1-ц2)

и совокупностью следующих граничных условий в точках х = О ,

с УФ) = ■ Л^Й = КНЬ)) = = (2)

При этом функции g , И и / имеют вид

(»ш"4 > Акад=IX м*))"" > /и=I • (з)

В (1) О - изгибная жёсткость пластины; N - сжимающее (растягивающее) усилие; Г7, р0, а -

скорость газа, плотность и скорость звука, соответствующие однородному потоку; М - число Маха, а. (у -1 +°о) - коэффициенты, характеризующие жёсткость основания; интегральный член

учитывает нелинейное воздействие продольного усилия; соу' - член, учитывающий аэродинамическое воздействие; а0 =1 (а0 =2) соответствует одностороннему (двустороннему)

обтеканию пластины; ы(х) - прогиб пластины; Е - модуль упругости; и - коэффициент Пуассона; Т7 - площадь поперечного сечения: J - момент инерции сечения. Все коэффициенты, входящие в уравнение и граничные условия, постоянные. В (2) коэффициенты сп й} (г = 0^л3у = 0 + т) -

произвольные, часть из них может равняться пулю; в зависимости от значений этих коэффициентов условия могут быть или линейными, или нелинейными. Значения т и п могут быть равными оо . Уравнение (1) возьмем в виде

/)м/4) + аи>' + аУ -Ы' \(м>')2 ёх = 0 . (4)

о

Перейдём в уравнении (4) к безразмерным переменным. После замены х - 1х , = , где I -некоторый характерный размер, а величины с чертой - безразмерные переменные, приходим к уравнению

^<4>+——= 0. • (5)

Будем изучать решения уравнения (5) при следующих граничных условиях:

н7"(0) - 0 , Й7'"(0) = 0Щ) = , п'О) = ^"(1) • (6)

Эти условия соответствуют свободному концу слева и упругой связи на правом конце пластины. В дальнейшем черту над переменными опускаем. Задача (5). (6) решалась численно. Численная реализация заключается в сведении краевой задачи к начальной задаче Коши. Уравнение (5) четвёртого порядка, а в исходной постановке имеются два начальных условия, поэтому для начальной задачи Коши не хватает двух условий. Запишем эти условия в виде м'(0) - X , м/(0) -V ,

П. А. Вельмисов, С. В. Киреев, 2004

где А,, V - параметры. Тогда м-(1) и м/(1) являются функциями X и V : F, (X,\') = w(l)-w"'([), F2(X>v) = = м>'(1) - м/"(1). Решаем следующую задачу Кош и:

а?2, .

ч>{л) +

аЛь з е^3

в

И' -

В

|(м/)2 с1х = О,

о

(7)

ъ(0) = Х, м/(0)=У, м/'(0) = 0, м/"(0) = 0

Задача Кош и (7) будет соответствовать краевой задаче (5), (6). если выполнятся условия = 0 и (Я,,у) = 0. Параметры А, и V будем определять с помощью ньютоновского процесса

по формулам

/л \

А

у

V

Г Ж (Ч .V „) др\ (Х„ .V „)

Л"1

дХ

О

ОУ

V

Ч п

^ (К »V „ )

) Ъ¥г (Хп ,у „) ЙК, (А.,, ,у „) {F2(X,l,vJJ

(3)

У

<8 И

ЭХ ЗУ

Этот итерационный процесс будем продолжать до тех пор, пока не выполнятся условия Р \Е2\<г у где 8 - заданная точность вычисления. Введём следующие обозначения. Пусть =

= и/' = у3 , VI/" = уА9 тогда интегро-дифференциальное уравнение в (7) можно записать в виде системы

1 1

у[ = У 2, У; у2 - /п ■ у! + ы* / • , где I = \tWfdx = \y\dx. (9)

Обозначим У -

У 2 у = 5 0 V

У> 0

КУ А у

\

, Г) =

г

Уъ

У4

(Ю)

ал6 , е^3

л

I)

за

у

Тогда задача Коши (7) примет вид

{У - ^(х, 7)

(П)

Г(0) = П

Задач}; Коши (11) решаем методом Рунге-Кутта шестого порядка с контролем погрешности на шаге. Сложность задачи Коши заключается в том, что в уравнении (11) присутствует интегральное слагаемое, для вычисления которого требуются значения подынтегральной функции сразу на всем отрезке интегрирования, что делает невозможным прямое применение метода Рунге-Кутта. Поэтому значение интеграла в уравнении (11) может быть получено из следующего итерационного процесса. Представим решаемое уравнение в виде

ы-

%) + %) = 0,где = и2ск, ) - и'

Г) о

(4)

, аъг з

—-щ ' к==1>2>-

1. Решаем уравнение О^^^О методом Рунге-Кутта по формулам

кл = к((х,у), /с2 = /?/

/

V

/2 /С

Х + -, у + — 2 2

\

/

/

ь

1

\

, къ =/?/ х + -5у + -(/с, +/с2) , кА =1г/(х + 11,у-к2+2къ),

2 4

V

/с5 = А/

/

х +

21г

\

1 ^ Г

27 у V

к

х + -,у+—(28/с, -\25кг -г54б/с3 + 54/с4 -378£5) 5 ' 625

Л

у

АУ = -(/с1 + 4/с3 + А',). 7(х + К) - 7(А) = г + 0(к6), г = — (42*, + 224А:3 + 2\кл-\62/с5 - 125/сб). (12) 6 ^ 336

Получаем значения и/ (х) и (к -1) на всем отрезке интегрирования, необходимые для

вычисления интегрального слагаемого /(ж,);

1.90728213151345Е+0000 1.78738В55905567Е«-0000 1.Б6549434839083ЕЧЮ00

-1.В654В^3^839083Е^0000 ■1.787386559055Б7Е-0000 ■1 90728213151345Е+00СС

1872 428 1868.121 1663.825

-017 ■0.172

-0.174 -0175 ■0.174 40) -0.11

--01«

-0.114

-0,114 -01« -0.1»

♦ (0) * - С ДОМ»? 13770182

♦ (1) « - 0 1737М8979541&

- 1 0 19 0 1:.; 0.1« 0 1 м 1 —

- - 01ГЗ ' . ^ II МВ

[ и 5.1 с -0 0.17с 01?4 0173 017 - 1

03 с;

1(0) = 0.119Л37701 «2 1(1) « 0.17'У7М&97!г541Ы

Рис. 1. Бифуркационные диаграммы

Рис. 2. Прогиб пластины для асимптотического решения

£ ~

= -0.01

0.458

-1 73930094789766Е-0001

-1.82859173048468Е-0001 -1.84014794706152Е-0001

1.84014794706152Е-0001 1.82859173048468Е-0001

1.79930094789766Е-0001

Рис. 3. Прогиб пластины при численном решении

2. Находим 1(м>{). Выражение /(м>,) содержит интеграл, который вычисляем с помощью

ь

квадратурной формулы Ньютона-Котеса. Пусть вычисляется интеграл ^(х)с1х. Отрезок

а

интегрирования \a\b\ разделим на «одинаковых частей длины к = {Ь-а)!п. Число п выбираем

кратным 5, для того чтобы весь интервал интегрирования разбился на участки, на которых подынтегральную функцию д(х) будем аппроксимировать интерполяционным многочленом

5 5

Лагранжа 1$(х) четвёртой степени д(х)«15 (х) = ^ д(х1)р. (х), где р. (х) -

весовая

/=0

;=0 Х1 Х; 14/

ь

функция. Тогда интеграл ^(х)с!х будет определяться в виде суммы интегралов вида

а

*5 5

х»=5к{р0д(х0) + рд{хх) + р2д(х2) + ръд{хъ) + Ра9Ы + Р5д(х5)), (13) 19 75 50

3. Решаем уравнение /)(ж,) + /(н-'.) = 0 методом Рунге-Кутта (12). Получаем значения ы'2(х) и Ц'(х) (к = 2) на всем отрезке интегрирования.

4. Находим 1{м>2), используя формулу (13).

5. И т. д.

Итерационный процесс продолжаем до тех пор, пока не выполнится условие тлах|>^ (х) - wk_< (,x)j < е , где 8 - то же самое, что и в ньютоновском процессе; wk (х) - прогиб на к -м

шаге; (х) - прогиб на (к - 1)-м шаге. Программа написана на языке Delphi 6. С помощью этой

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

л/м7 — 1 - Via => а -а0р0aV . Если возмущение е = X - si образуется за счёт увеличения скорости потока, то ветвление в. точке надкритическое. Если 8 образуется вследствие уменьшения

• » i • I *

• ... . У \ V • - Ч

изгибной жёсткости, то ветвление в точке Х0 =s}0 будет подкритическим. На рисунке 1 представлены

бифуркационные диаграммы прогиба пластины при фиксированных коэффициентах изгибной жёсткости А < Д < Д в зависимости от изменения скорости потока сверх критических значений

Х} < Х2, =sl < Х^. На рисунке 1 крайний левый график соответствует \ , Д, средний - X,, Д, крайний правый - ХЪ9 Д . Графики построены при: 1, р0 = 330, 0 =1 Н/м, I -10 м,

Д = 106,999 кг/м2, Д = 10' кг/м2, Д =\01ОО] кг/м2. На рисунке 1 верхние ветви диаграмм

соответствуют положительным решениям, а нижние - отрицательным. Была проведена проверка полученных численных решений с аналитическими с использованием математического пакета Matbcad 2000 Professional. С помощью этого пакета вычислялись коэффициенты, входящие в асимптотическое решение, полученное аналитически методом Ляпунова-Шмидта. На рисунке 2 представлены формы прогиба пластины асимптотического решения, где ф(х) соответствует

положительному решению (ф(*) = w(x)); и/ (х) - отрицательному (ф(*) = —vi>(*)), а на рисунке 3 ~

численного решения. Точность совпадения асимптотического решения в первом приближении с

численным решением составляет порядок 10" .

• •••••»••••••••••••с аэрогидромеханике, аэрогидроупругости,

математическому моделированию. Велъмисов Петр Александрович, доктор Киреев Сергей Владимирович, ассистент

физико-математических наук, профессор, кафедры «Высшая математика » УлГТУ. Имеет заведующий кафедрой «Высшая математика» статьи по аэрогидроупругости,

УлГТУ. Имеет монографии и статьи по математическому моделированию.