Конечноэлементная математическая модель динамики криволинейного трубопровода с пульсирующим потоком рабочей жидкости Текст научной статьи по специальности «Физика»

УДК 533.951

КОНЕЧНОЭЛЕМЕНТНАЯ МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ДИНАМИКИ КРИВОЛИНЕЙНОГО ТРУБОПРОВОДА С ПУЛЬСИРУЮЩИМ ПОТОКОМ РАБОЧЕЙ ЖИДКОСТИ

© 2009 Т.Б. Миронова

Самарский государственный аэрокосмический университет,

Поступила в редакцию 27.04.2009

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

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

Ранее были разработаны математические модели процессов связанных колебаний жидкой среды и твердотельных элементов в элементарном объеме на границе раздела сред [3]. Однако данные модели слишком громоздки для решения задачи анализа виброакустических процессов в трубопроводах. Известна конечноэлементная модель динамики трубопроводов сложной пространственной конфигурации с пульсирующими потоком рабочей жидкости, созданная на базе использования программного комплекса ANSYS [1-2]. Недостаток данной модели заключается в ее высокой трудоёмкости. Указанный недостаток данной модели устранен в математической модели [3] Существенное сокращение вычислительной трудоемкости достигается при рассмотрении трубопровода с точки зрения механики стержней [4].

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

В настоящее время многие известные универсальные CAE-системы базируются на мето-

Миронова Татьяна Борисовна, аспирант кафедры "Автоматические системы энергетических установок". E-mail: mirtb@rambler.ru.

де конечных элементов (МКЭ), т.к. метод конечных разностей (МКР) считается недостаточно эффективным и в значительной мере устаревшим. МКЭ по сравнению с МКР требует меньше машинных ресурсов (меньше оперативной памяти), расчет идет быстрее (меньше затраты процессорного времени), результат расчетов более адекватен. Преимуществом конечноэлемен-тной модели перед конечно-разностной является также возможность внедрения разрабатываемых конечных элементов в современные универсальные САЕ-систем8, такие как Ansys, Nastran, Patran, Cosmos и др. что позволяет расширить их элементную базу, а также в частности решать задачи моделирования виброакустических характеристик трубопроводных систем с пульсирующим потоком рабочей жидкости в комплексе с анализом динамической нагружен-ности присоединенных агрегатов и систем, анализом технического объекта в целом.

Система уравнений, описывающая малые колебания пространственно криволинейных трубопроводов с осевой линией, лежащей в одной плоскости, при их силовом нагружении пульсирующим потоком рабочей жидкости и кинематическом возбуждении со стороны присоединенных агрегатов и систем записывается в виде [3]:

1 ди,

Жз д е '

2 д 2и, д4и, д4и,

Жз—2---- nw- 1

дт де дт

д6и,

де дт

~nwXi

2 д2и,

дедт

де6

-(

p + nw +

Ж

-((p + nwW )жз +Жз ))

д 4U,

де4

2 3w 3w

Y = nw%3--nw—,

де дт де

1 д и, Ф3=— -г-г + Жзи1 Ж де

где В - безразмерная координата, отсчитываемая вдоль линии центров тяжести сечения трубопровода от начала отсчёта до некоторого произвольного поперечного сечения; т - безразмерное вре-

мя; n = ■

m1 (s ) - погонная масса трубо-

провода; m2 (s) - погонная масса рабочей жидкости в трубопроводе; c3 - кривизна осевой линии в плоскости, перпендикулярной e3 ; e3 (s,t) - единичный вектор, направленный по бинормали к осевой линии трубопровода; e1 (s,t) - единичный вектор, направленный по касательной к осевой линии трубопровода; ё2 (s,t) - единичный вектор, направленный по нормали к осевой линии трубопровода; l - длина трубопровода; w - вектор безразмерной скорости рабочей жидкости; p - безразмерное давление; u1 - виброперемещение в направлении e1 (s,t); u2 - виброперемещение в направлении e1 (s, t).

Выражения для мгновенных значений колебательных составляющих давления и скорости рабочей жидкости в рассматриваемом случае гармонических колебаний в безразмерных величинах давления и скорости рабочей жидкости, представлены в виде (в левой части приведенных ниже равенств - безразмерные величины, а в правой - размерные, кроме t):

2Ф(-в) 2 ___

COS- 17/2

pM = - C Рвх 22

COs

2ф

B^coS 2ф2\т^,

w{e,r) = -

in2n(1 ~в)

2nfl

MJmi+msiJi2nfi2 l'"'+m2'

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

Рассмотрим схему решения второго уравнения системы (1). Данное уравнение является линейным нестационарным дифференциальными уравнением. Его можно записать в общем виде [5]:

du

d2u,

Lui + p-a—1 -P—t = 0 , на Q , (2)

dr

dr2

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

Для решения данного уравнения применим метод частичной дискретизации, при котором исходное дифференциальное уравнение с част-

ными производными заменяется системой обыкновенных дифференциальных уравнений [5]. Полученная система решается повторной дискретизацией по времени.

Аппроксимация и1 для решения и1 с помощью метода частичной дискретизации записывается в виде:

ui * ui = Z am(*)Nm(e) ,

(3)

где Нт(в) - базисные функции. При этом N т (в) не обязательно удовлетворяет всем краевым условиям.

Аппроксимирующее уравнение по методу взвешенных невязок в общем виде записывается:

J WlRQdQ+J WlRrdr =0,

Q Г

(4)

du1 0d2u1

где Ка = Ьи1 + р —а—---р , 2 - невязка апп-

dт dт

роксимации по области, ЯГ = Ми1 + г - невязка аппроксимации в краевых условиях, -линейно независимые весовые функции.

С учетом приведенных соотношений перепишем второе уравнение системы (1) в виде:

С д6й, ( 2 „ д4й,

J (-двй1 + ( + + Х — +

(( Л 2 4 \д 2щ + ( + nw — + Хо )—Г - nwXo — +

dw

де2

дт

dw 2 д й, д й, + nw--х30-2- + —Тг +

де дт де дт

(5)

+ nw

д4йх 2 д 1й1

де дт

+ nwX30

дедт

)Wlde = 0

Используем аппроксимацию по Галеркину. В этом случае вместо привлечения новой системы функций в качестве весовых множителей выбираются сами базисные функции:

= Nl и = - Nl на Г .

Подставим аппроксимацию и1 в уравнение (5) и запишем полученную систему дифференциальных уравнений в векторной форме:

м t -u]+c m+ь f ]=[f ];

(6)

Здесь компоненты отдельных матриц и правой части определяются равенствами:

м,

C,.

■ = J (""¡Т -Ч-N-)

t (d3Nm 2 dNmV..

=Jnw[-e+z*O-NTJN-e

Km =j

-6Nm ,„ 2 2-4Nm

—твт + (2Xo+P+nwf)—rir + de de

m

2

mt + m2

m=0

33

33

c

c

+(х30+р+

(е2

, г 2 (д™ д™„7 ,

* = 1пХ [-¡¡Г- ™ -дв] м°е

I--I'

^дг де)

В общем случае аппроксимация й1 на элементе с р +1 узлами будет сводится к многочлену степени р . На таком элементе узлы которого с номерами от 0 до р помещены в точки е0,

В1 , В2

е

ер, ассоциируемая с узлом I.

Базисная функция элемента Nl будет многочленом степени р , принимающим значения нуль во всех других узлах элемента.

В уравнение (5) входит производная 6 порядка. Следовательно, для получения точного решения базисные функции, входящие в аппроксимацию должны иметь 6 порядок или выше. Это необходимо для выполнения требования полноты системы базисных функций, позволяющей им с любой степенью точности аппроксимировать неизвестную функцию [5]. В качестве базисных функций удовлетворяющих данным условиям возьмем многочлен Лагранжа степени р = 6 [5]

Щ =л =\е-ео Хе-еМе-В-1Хе-е+1 Х—-ер )]х х[(е -в \е-е1 Х)е ^Лв еМе -ерМ (8)

Запишем базисные функции для одного се-миузлового элемента. Выразим их через нормированную локальную координату V , определяемую равенством:

2(е-еес)

V-

ке

(9)

Здесь еес - координата центра элемента, ке - длина элемента. Элемент принадлежит отрезку -1 <%< 1:

N1 = —(V6-V5 - -V + -V + —V2 - —V);

0 80 9 9 81 81 '

N = --2V--10V—+20% +1V -2;

1 40 3 9 27 9 27 '

N = —V6 -1V5 -13V4 +13V3 +12V2 - — V);

2 16 3 9 27 27 27 '

N =--(V6--V4 + --);

3 4Ь 9* 81 81'

(10)

N4, = — (V6 +-? -—V -—V +-¥ +—?);

4 16 * 3* 9* 27* 9* 27 ЛТ, 243 6 2 .5 10 .4 20 ,3 1 ,2 2

N. =--(V + —V--V--V +-£2 +—V);

5 40 3 9 27 9 27

N1 = —(V6 + {5 - -V4 - -V3 + — V2 + — V).

6 80 * 9* 9* 81 81

Полученные базисные функции представлены на рис. 1.

Запишем систему уравнений (7) через нормируемую локальную координату V . Получим:

м„

с„

-I

2 d2Nm 2 к'

- х2 —N к' (IV *30 2 г

N Ое

1 ((2 ^2dзNm

=НН +х

2 dNm

N¡dе

К1т =|

2 Л5 о6 N

2 Л3 Сы

+ (Х0 + Р + п™2)Х 1-2-1

(V3 30 (V +(2Х+р) "(—г + (11)

тс ^

^ к') О

, 1 2 (д™ д™) к'т1е /,=1пХ (дт- ™ д—) т™

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

Рис. 1. Одномерный элемент и ассоциируемые с ними стандартные базисные функции шестой степени

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

а « а = 2^т(т) •

(12)

Базисные функции Nт (т) должны иметь степень не ниже второй, так как в уравнение (6) входят вторые производные по времени.

Возьмем типичный квадратичный элемент п по времени с тремя узлами, помещенный в точки т2п , Т2п+1 , т2п+2 (см. рис. 2 3).

На этом элементе п имеем

а = а2п^п + а2п+1^п+1 + а2п+2Щп+2. (13)

Тогда как все остальные базисные функции на элементе п равны нулю. Здесь

+[-2М+(1-2ГЖЯДО +(1/2-2^+^X40]] +(16) +\м-(1-№пС(Ю+(1/ 2+.р-г^жк)]" = /Ж2

Интерполируя / тем же самым способом, что и а , получим

/п =&2п+2 + (1/2 - 2р+у)/2п+1 + (1/2+Р-у)/2п .(17)

Рассмотрим частный случай краевых условий - жесткую заделку концов трубопровода. В этом случае при в = 0 и в = 1 и1 = и2 = 0 . Тогда из первого уравнения системы (1) следует, что

ди

1 = 0 ТОТ™„ д и1

, а из третьего уравнения

= 0. Опу-

дв ' ^ —...... дв2

стим индекс 1 у параметра и.

Для рассматриваемой конечноэлементной модели граничные условия можно записать в виде

Nn =

й2^ dт2

- Т(1 - Т) -1/2 + Т

2

2 п

йт

А

А

; ¡V" =1-Т

' 1 у 2п+1 1 1

2

йт

2Т А ;

___^ = Т(1+Т) ймпп+2 = 1/2+Т

А2 ; *2п+2 = 2

йт2

йт

А

й2^пп+2 = _Т = г - г2п+1 йт2 А2 ; Ап

А =т, +,-т, .

2 + 1 2

; Ап =т 2 п+2 т 2п+1 ;

(14)

2 ит = 0

т=0

М-1 ди

т=0

дв

= 0 при в = 0 и в = 1.

(18)

М-1 д 2и 2—т=0 ^ дВ2

т=0

Распишем последние два уравнения системы (18), подставив в них аппроксимацию (3) и учитывая, что перемещение в узлах на границе а0 и а6, а также базисные функции N1, N2 , N3 , N4 , N,5 в граничных точках равны нулю а0 = 0 и а6 = 0 , получим

Применение к уравнению второго порядка (6) стандартного метода взвешенных невязок дает

М—- + С— + Ка - / йт йт

Wndт = 0. (15)

п = 0,1,2,...

Здесь интегрирование производится только по элементу , поэтому в полученном выражении можно подставить для а значение (13). Учтем, что матрицы С и К зависят от времени. Тогда в силу (14) уравнение метода взвешенных невязок после соответствующих преобразований принимает вид

Рис. 2. Разбиение временной области на квадратичные конечные элементы

дЖ(0) д#2(0) д#3(0) а —^^ + а2—+ а3—+ дв дв дв

+ а

дNA{0) дNъ{0)

4 - + а5

дв дв

= 0

дЖ(1) д#2(1) дЖ(1) а —+ а2 —+ а3—+

1 ^ 2 ^ 3

дв дв дв

+ а

авд , д#5(1)

4 " + а5 дв дв

= 0

д2 N,(0) д2 #2(0) д2 #3(0)

а1-^ + а2-+ а3-^ +

дв дв дв

+ а

I2 N4(0) , а д2N5(0) = 0

дв

+ а

4 2 5 2

дв2

а

д2ад), ^ д2N2(1), „ д2N3(1)

- + а

+ а

1 дв2 2 дв2 3 дв2

- +

+ а

д 2 N4®.+а, д:М)=0

4 де1 5 де1

т=1

г

т

Рис. 3. Базисные функции элемента п Таким образом имеем 4 уравнения и 5 неиз-

вестных величин а,, а,, а,, а,, а.

а3 нахо-

Ч ' "2' 3' **5

дим из уравнения (6), которое для случая одно го элемента можно записать в виде

[м„ ]+[с„ (

От

+ \К 33 [3 ]=[ ].

(20)

Необходимо отметить, что изменение граничных условий не приводит к перестройке ко-нечноэлементной схемы решения. Выражения (11) для внутренних узлов элемента остаются прежними. Изменяется только вид уравнений системы (19), описывающих граничные условия. В конечно-разностной модели изменение граничных условий приводит к перестройке схемы решения, что увеличивает трудоемкость получения результатов расчета по сравнению с предложенной конечноэлементной моделью.

В качестве модельного примера возьмем трубопровод со следующими параметрами: I = 0,4м ; О = 0,004м ; 8 = 0,0006м ; р = 7800кг / м3, Е = 2 • 1011 Па; Я = 0,23 м ; где О - наружный диаметр трубопровода, 8 - толщина стенки трубопровода; Я - радиус кривизны.

Предположим, что в начальный момент времени деформация трубопровода отсутствует и он находился в состоянии покоя. Затем трубопровод нагружается установившимися колебаниями рабочей жидкости в которой реализуется стоячая волна с параметрами / = 250Гц, рвх = 105Па. Характеристики рабочей жидкости : рж = 870кг/м3, с = 1300м/с.

При проведении расчётов возьмем один се-миузловой элемент и постоянный шаг дискрети-

зации по времени М = —. В качестве начальных 70

данных в рассматриваемой трехслойной схеме для начала вычислений требуются начальные данные, которые задаются в виде

Оа

а(т = 0) = а0 = 0 , ~^(т = 0) = 0

(21)

Для определения а] применим стартовую ко-

Оа а'+ - а'

нечно-разностную схему — = -

М

-, где ' - со-

ответствующий временной слой.

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

в> (1/2+г)2 , у> 2, 2 + Г + в> 0. (22)

Значения в и У соответствуют весовым функциям различных видов. Были выбраны следующие значения: в = ~ и у = 1. При таких значениях схема является безусловно устойчивой и не дает искусственного числового затухания.

Результаты расчета представлены на рис. 4, 5.

При анализе рис. 4 обращает внимание негармонический характер вибрации при возмущении гармонической волной давления. Это связано с тем, что предложенная модель построена без учёта сил трения. Свободное движение в расчётной схеме не затухает, накладываясь на вынужденное движение системы, что и обуславливает негармонический характер вибрации на расчётных временных реализациях. Сказанное подтверждает рис. 5 на котором показаны амплитудные спектры, соответствующие рис. 4. На амплитудных спектрах явно видны два пика -250 и 341 Гц, первый из которых соответствует вынужденным колебаниям, а второй - собственным колебаниям по первой собственной форме.

Собственные частоты, полученные по разработанной конечно-элементной модели, сравнивались с собственными частотами, полученны-

Рис. 4. Расчетные временные реализации нормального вибросмещения при колебаниях трубопровода:

а - в 1 узле е = 0,33; б - в 3 узле е = 0,5

Рис. 5. Амплитудные спектры нормального вибросмещения при колебаниях трубопровода:

а - в 1 узле е = 0,33; б - в 3 узле е = 0,5

ми в программном комплексе А^УЗ и полученными с помощью конечно-разностной модели. Для рассматриваемого модельного трубопровода, собственная частота полученная в А^УЗ равна 336 Гц, по конечно-разностной модели 340 Гц, по конечноэлементной модели - 341 Гц.

Система уравнений (1) получена при неучёте сил трения (сопротивления) в материале трубопровода и его опорах. Силы трения в опорах задаются соответствующей формулировкой краевых условий. Введём в первое уравнение системы (1) член, характеризующий действие сил трения (сопротивления).

В работе [3] этот член записывался в виде

- Н

дй1 ИТ

. Данное выражение описывает сопро-

тивление малых свободных колебаний линейной системы, обобщенная сила которого пропорци-

дй, „

ональна скорости —-, и Н > 0. Эта модель тре-дт

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

и равна Н

д 5й1 де4дт

С учётом этого второе уравнение системы (1) можно переписать в виде:

2 д2щ

1за

д4и1 д4м1

дт

2

-нх д и

де2дт2

- пм-

деъдт

2 д\ Мзо-1

Сзо

дТде деь

-+

дедт

д4и,

((+пМ)+2Хо)) + (23)

+

д2и1

(р+пМ)

,2^2, 2 дМ

+Хзо Дзо +

де

де

2 дм

Х-Г. дт

Моделирование виброакустических характеристик с учетом сил трения проведем для ранее рассмотренного трубопровода при тех же условиях его нагружения.

На рис. 6 в безразмерных параметрах изображены расчётные временные реализации нормального вибросмещения в узловых точках. Рис. 6, а соответствует системе без трения (н = 0), в которой наряду с установившимися вынужденными колебаниями реализуются незатухающие свободные колебания. Введение в систему трения (н > 0) приводит к затуханию свободных колебаний (рис. 6, б, в). С течением времени расчётная временная реализация вибросмещения переходит в синусоиду с частотой пульсаций давления. При этом, чем больше Н, тем быст-

Л 1

Л \И м г

1 I и V у 1

t

а) б) в)

Рис. 6. Временные реализации безразмерного нормального вибросмещения в 1 узле е = 0,33 при вынужденных колебаниях трубопровода под действием пульсирующего потока жидкости из положения покоя при различных величинах параметра трения: а - Н=0; б - Н=0, 003; в - Н=0,03

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

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

На рис. 7 представлена амплитудная частотная характеристика виброускорения сечения трубопровода, колеблющегося под действием пульсирующего потока рабочей жидкости с амплитудой 105 Па во входном сечении, полученная в результате расчетов по разработанной модели. При расчетах параметр трения Н принимался равным 0,03. Обращает на себя внимание наличие резонансного пика в районе частоты £=341 Гц. Он связан с резонансом механической подсистемы трубопровода на этой частоте и соответствует первой собственной частоте данного трубопровода, полученной выше при моделировании его собственных колебаний. Второй резонансный пик на частоте 812,5 Гц соответствует четверть волновому гидравлическому резонансу.

В работе представлена конечноэлементная модель, описывающая динамические характеристики пространственно криволинейного трубо-

провода при его силовом нагружении пульсирующим потоком рабочей жидкости. Рассмотрен частный случай решаемой задачи - вибрация трубопровода, ось которого лежит в одной плоскости, под действием стоячей волны в рабочей жидкости. Сравнение результатов расчетов по предложенной модели с результатами, полученными с помощью конечно-разностной модели [3] и с результатами, полученными с помощью программного комплекса ANSYS, подтвердили адекватность модели. Кроме указанных выше, еще одним преимуществом разработанной модели по сравнению с конечно-разностной заключается в том, что при изменении граничных условий не требуется перестройка расчетной схемы. Это снижает трудоемкость получения решения по предложенной конечноэлементной модели по сравнению с конечно-разностной моделью.

СПИСОК ЛИТЕРАТУРЫ

1. Макарьянц Г.М., Прокофьев А.Б., Шахматов Е.В., Шестаков Г.В. Исследование виброакустических характеристик трубопровода при его силовом нагружении с использованием программного комплекса ANSYS // Сборник трудов четвертой конференции пользователей программного обеспечения CAD-FEM GmbH. М.: Полигон-пресс, 2004. С. 280-287.

2. Макарьянц Г.М. Разработка методик расчета и исследование виброакустических характеристик трубопроводных систем. Дисс... канд. техн. наук. Самара: СГАУ, 2004. 191 с.

3. Прокофьев А.Б. Разработка метода комплексного анализа динамики и прочности трубопроводных систем с гасителями колебаний рабочей жидкости. Дисс... докт. техн. наук. Самара: СГАУ, 2008. 191 с.

4. Светлицкий В.А. Механика стержней. М.: Высшая школа, 1987. 304 с.

5. Зенкевич О., Морган К. Конечные элементы и аппрок-симация/М.: Мир, 1986. 318 с.

6. Бабаков И.М. Теория колебаний. М.: Дрофа, 2004. 591 с.

FINITE ELEMENT MATHEMATICAL MODEL OF PIPELINE DYNAMICS UNDER FORCE EXCITATION BY A PULSATING FLOW

© 2009 T.B. Mironova

Samara State Aerospace University

The finite element mathematical model of vibroacoustical characteristics of pipe system is developed. The model allows calculations of pipe characteristics under force excitation by oscillating fluid flow, including unsteady conditions. The pipe system can be of complex configuration. Results of calculation on the developed model are presented.

Keywords: pipeline system, dynamic characteristic, power fluid, vibration, finite element model, fluctuation.

Tatiana Mironova, Graduate Student at the Automatic Systems of Power Plant Department. E-mail: mirtb@rambler.ru.