CC BY
10
ЕСТЕСТВЕННЫЕ НАУКИ
УДК 539.3:533.5:517.9
П. А. ВЕЛЬМИСОВ, С. В. КИРЕЕВ
ЧИСЛЕННЫЙ МЕТОД РЕШЕНИЯ ЗАДАЧИ О СТАТИЧЕСКОЙ НЕУСТОЙЧИВОСТИ ТРУБОПРОВОДА
Предложен численный метод решения нелинейной краевой задачи о статической неустойчивости трубопровода с протекающей в нём жидкостью. Получены бифуркационные диаграммы, показывающие зависимость максимального прогиба трубопровода от скорости протекания жидкости.
Математическая модель задачи об изгибных формах трубопровода, с протекающей в нём жидкостью, описывается нелинейным интегро-
дифференциальным уравнением
+
и™" + /(^ - е»'" |(и/)2 йх -о,
о
ЕР
N = тМ , 0 = — ,£> = £/
2£
(О
и совокупностью следующих граничных условий в точках х - 0 , х = Ч
с уф)=, ау\ь)=н{ут,
Ь = 0,Ь = £.
При этом функции g , к и / имеют вид
(2)
п
к=\
т
А(н< ,
X
2/)+1
(3)
//-О
В (I) О - изгибная жёсткость трубопровода; >0
- сжимающее (< 0 - растягивающее) усилие; т*
- удельная масса жидкости; N - сжимающее (растягивающее) усилие; II - скорость движения жидкости;
а. О = 1 ~ коэффициенты, характеризующие
жёсткость основания; интегральный член учитывает нелинейное воздействие продольного усилия; ы(х) -
прогиб трубопровода; Е - модуль упругости; Р -площадь поперечного сечения; У - момент инерции сечения. Все коэффициенты, входящие в уравнение, и граничные условия постоянные. В (2) коэффициенты
П. А. Вельмисов, С. В. Киреев, 2005
с:, с1. (/ = 0 -г- п, у = 0 -г- т) - произвольные, часть
из
них может равняться нулю; в зависимости от значений этих коэффициентов условия могут быть или линейными, или нелинейными. Значения т и п могут быть равными оо .
Уравнение (1) возьмём в виде
. М?2 _„. а/ _з
Ху<4)+—й/3-
о
о
= О,
о
м_
~Ъ
(4)
где - некоторый характерный размер, а величины с чертой - безразмерные переменные (х = £х ,
Будем изучать решения уравнения (4) при следующих граничных условиях: \¥'(0) = 0,
= 0,
м7(1) = 0, Щ\) = 0.
(5)
Эти условия соответствуют свободному элементу слева и жёсткой заделке на правом конце трубопровода. В дальнейшем черту над переменными опускаем. Задача (4), (5) решалась численно. Численная реализация заключается в сведении краевой задачи к начальной задаче Коши. Уравнение (4) четвёртого порядка, а в исходной постановке имеется два начальных условия, поэтому для начальной задачи Коши не хватает двух условий. Запишем эти условия в виде и^(0) = X ,
\н>"(0)=у, где А., V - параметры. Тогда Ц1) и
м/(1)
являются
функциями
к
и
V :Р{(Ху)-=ЕИ<1Д,у), ^лО^ХилО. Решаем следующую задачу Коши:
(4) №2 з
И' +-Н' +-VI' -
О
й
ег о
I
О
= ^ и/(0) = 0, У/(0)=У, = (6)
Задача Коши (6) будет соответствовать краевой задаче (4), (5), если выполнятся условия Г](к>у) = 0 и F2(X,v) = 0. Параметры ¡1 и V будем определять с
помощью ньютоновского процесса по формулам
л \ и
V 'Н-1У V ,
а
-1
дк
Этот итерационный процесс будем продолжать до тех пор, пока не выполнятся условия и /\ <е ,
где е - заданная точность вычисления. Введём следующие обозначения. Пусть ту = у{, = у,,
IV" = уъ, и " = уА, тогда интегро-дифференциальное уравнение в (7) можно записать ь виде системы
- У2=Уз> У',=У4>
у'А Ю-уъ/П-у1 у319
I I
где /=
(В)
о
о
Обозначим
V
>2 >^0 = 0 V
V«,
Г
Пх, У) =
У 2 Уг У л
\
N1
Уз
з
ю3 г
К У г1
(9)
\
У
(10)
Тогда задача Коши (6) примет вид
{У = Р(х, У) У(0) = У0.
Задачу Коши (10) решаем методом Рунге-Кутта шестого порядка с контролем погрешности на шаге. Сложность задачи Коши заключается в том, что в уравнении (10) присутствует интегральное слагаемое, для вычисления которого требуются значения подынтегральной функции сразу на всем отрезке интегрирования, что делает невозможным прямое применение
метода Рунге-Кутта. Поэтому значение интеграла в уравнении (10) может быть получено из следующего итерационного процесса. Представим решаемое уравнение в виде
е*3
/>Ю + /( = где / = -
О
к
|и\с1х,
о
МР 2 п Р6
) = + V, к = 1,2,...
о л й
1. Решаем уравнение = 0 методом Рунге-
Кутта по формулам
к к 4
(7) к=к/(х,у),
к2 = к/
\
х + -,у + — 2 2
/
/с, + + + £,)
= к/(х + Л, у - к2 + 2к}),
/
= Ы
/
л
v
27
/
¿6 - Л/
/ /2
V 1
5 625 1-34
5
6
+ /г)-У(/2) = г + 0(Л6),
1
г =
336
(42*, +224£3 +21£4 -162Д:5 -125/с6). (11)
Получаем значения и уу"(х) (& = 1) на всем
отрезке интегрирования, необходимые для вычисления интегрального слагаемого ).
2. Находим /(VI/,). Выражение /(>',) содержит интеграл, который вычисляем с помощью квадратурной
формулы Ньютона-Котеса. Пусть вычисляется инте-
ь
фал [п(х)с1х. Подынтегральную функцию д(х) бу-
и
дем аппроксимировать интерполяционным многочленом Лагранжа четвёртой степени
1
д(х)» ¿5 (*) = £ )/?. (*),
где
1=0
5 X - X :
4 Л-(*) = --~ ~~ весовая функция. Тогда инте-
М Х1 ~
о
грал будет определяться в виде суммы инте-
л
гралов вида
х5
х5
+
Jq(x)dx « Jls (x)dx =£Pi q(Xi) =5h(p0q(x0)
xo xo i=0
+p1q(x1)+p2q(x2)+p3q(x3)+p4q(x4)+p5q(x5)), (12)
19 75 50
где p0 = p$ = —, /7, = />4 = — , p2 = ръ =
max
288' * ' 288 " " 288
3. Решаем уравнение D(w2) + /(w,) = 0 методом
Рунге-Кутта (11). Получаем значения и/(л*) и н>''(я:) (к -2) на всём отрезке интегрирования.
4. Находим I(w2) используя формулу (12).
5. И т. д.
Итерационный процесс продолжаем до тех пор, пока не выполнится условие
wk (х) - wA4 (jc)| < S , где 8 - то же самое, что и в ньютоновском процессе; wk(x) - прогиб на /с-ом шаге; - прогиб на (&-1)-м шаге. Программа,
реализующая алгоритм решения задачи, написана на языке Delphi 7. С помощью этой программы строятся бифуркационные диаграммы, показывающие зависимость максимального прогиба трубопровода от скорости протекающей жидкости. Если возмущение
8 = Х* - si образуется за счёт увеличения скорости потока, то ветвление в точке Х0 = s* надкритическое.
Если 8 образуется вследствие уменьшения изгибной
жёсткости, то ветвление в точке Х0 = s^ будет под-критическим. Рассмотрена модель с параметрами (Е = 206 • 109Н/м2 (сталь), ¿ = 1м, 0=35-1О3Н/м, а3 =1, D2 = 449 Нм2, Do=450 HM2,
D3 = 451 Нм2, А^о = 1Н, от, =10кг).
Для неё построены бифуркационные диаграммы прогиба трубопровода, представленные на рисунке 1 (построенные с помощью написанной программы) и 2 (построенные в Mathcad 2001 i Professional) при фиксированных коэффициентах изгибной жёсткости D, < Д < Д в зависимости от изменения скорости
потока сверх критических значений Х{ <Х2 = Sg < Хъ. На рисунке 1 крайний левый график соответствует \ , Dx, средний - Х2, Д, крайний правый - Ху, Д . На рисунке 1 и 2 верхние ветви
диаграмм соответствуют положительным решениям, а нижние - отрицательным.
0.05798652 0.05597066 0.05387823
0i
-0.05387823 -0.05597066 ■0.05798652
21.390002 U
ч
21 035477 21.072072 21.048640
Рис.1
D!
0025 О
-D025 -О 05 •С073
а\
0015 0.0135 0.012 0.0105 0 009 Ф(х) CD075
- 0 006
0.0045 0003 0 0015 0
г fl
i
1 i i
1 21W 21 iß 21 12 21.56 2t 2 2124 2! 2 2132 21 36 21 Рис. 2 r - 001
Ч)
о
<21 мщ .ocsm&
vj(2l иЗСф.ООЙММР* *}(21 ЭЮХП) =0йЯ92КЯ
».Y21202103). О 0321*322
о
-0 0015
-о.ооз
-0.0045 -0 00Ö V(xh00075
--0.009
-0.0105 -0.012 -0 0135 -0.015
«q> - 0.010208601 КО-о
v(0) - -O.OIG20S601
4*0-0
Рис.3
£•=0.01
w
0.01022629
0.01
I MBB
£=0.01
w
0.01022629
001022S2S
W
Рис.4
Была проведена проверка полученных численных решений с аналитическими с использованием математического пакета Mathcad 200Ii Professional. С помощью этого пакета вычислялись коэффициенты, входящие в асимптотическое решение, полученное аналитически методом Ляпунова-Шмидта. На рисунке 3 (построенные в Mathcad 200Ii Professional) и рисунке 4 (построенные с помощью написанной программы) представлены формы прогиба трубопровода асимптотического решения, где ф(лс) соответствует положительному решению (ф(л-) = w(x)), V|/(x) - отрицательному (ф(х) = ->ф0),при D0 .Ьо-
Точность совпадения асимптотического решения в первом приближении с численным решением составляет 99.83%.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Вельмисов, П. А. Численный метод решения задачи о статической неустойчивости пластины в сверхзвуковом потоке газа / П. А. Вельмисов, С. В. Киреев // Труды СВМО. - Саранск. - 2004. - Т. 6, № I. -С. 166-170.
2. Вельмисов, П. А. О статической неустойчивости трубопровода / П. А. Вельмисов, С. В. Киреев //Численные и аналитические методы расчёта конструкций: Труды международной конференции. - Самара.-1998.-С. 244-249.
Вельмисов Петр Александрович, доктор физико-математических наук, профессор, заведующий кафедрой «Высшая математика» УлГТУ. Имеет монографии и статьи по аэрогидромеханике, аэрогидро-у пру гости, математическому моделированию.
Киреев Сергей Владимирович, старший преподаватель кафедры «Высшая математика» УлГТУ. Имеет статьи по аэрогидроупругости, математическому моделированию.
\
