CC BY
12
УДК 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" .
• •••••»••••••••••••с аэрогидромеханике, аэрогидроупругости,
математическому моделированию. Велъмисов Петр Александрович, доктор Киреев Сергей Владимирович, ассистент
физико-математических наук, профессор, кафедры «Высшая математика » УлГТУ. Имеет заведующий кафедрой «Высшая математика» статьи по аэрогидроупругости,
УлГТУ. Имеет монографии и статьи по математическому моделированию.
