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

УДК 517.9; 539.3; 519.6

Численное моделирование движения абсолютно гибкого стержня в потоке воздуха

© Ф.Д. Сорокин, Ф.Р. Низаметдинов МГТУ им. Н.Э. Баумана, Москва, 105005, Россия

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

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

Введение. При проектировании и расчете электропередач, канатных дорог, гибких трубопроводов и других конструкций, содержащих провода или тросы, возникает проблема взаимодействия абсолютно гибкого стержня с внешним потоком жидкости или газа (далее — «потоком»). Поведение таких конструкций в случае нерастяжимых стержней описывает нелинейная система дифференциальных уравнений в частных производных [1]:

где р — плотность материала стержня; А — площадь поперечного сечения; г — радиус-вектор оси стержня; q — распределенная нагрузка со стороны потока и другие виды нагрузок; I — время; ^ — осевая координата.

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

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

При проектировании конструкций, взаимодействующих с внешним потоком, как правило, рассматривают следующие основные задачи [4-10]:

1) определение равновесной формы гибкого стержня в стационарном потоке;

2) исследование динамической устойчивости конструкции в потоке;

3) определение спектра собственных частот для отстройки от резонансных режимов при возможном кинематическом или силовом возбуждении.

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

Описание алгоритма. Для перехода к СОДУ стержень необходимо разбить на прямолинейные конечные элементы (КЭ), соединенные шарнирно в узлах конструкции. Участки стержня будем считать линейно-упругими и невесомыми, а их массы — сосредоточенными в узловых точках (рис. 1).

Рис. 1. Дискретная модель абсолютно гибкого стержня (цифрами показаны узловые точки)

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

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

стема — изменяемая или неизменяемая, статистически определимая или неопределимая. Достаточно, чтобы расчетную схему объекта можно было свести к набору невесомых стержней с сосредоточенными массами в узлах. Аналогичные модели рассмотрены в [11].

Рассматриваемый КЭ учитывает как упругие, так и инерционные свойства участка стержня (рис. 2).

ть

ЕА

'о ^^

Рис. 2. Конечный элемент абсолютного гибкого стержня Продольную силу в КЭ определяем в соответствии с законом Гука:

N = ЕЛ — = еЛ—^ , (1)

/ ¿0

где Е — модуль упругости материала стержня; Л — площадь поперечного сечения элемента; /0, / — длина КЭ в исходном и деформированном состоянии соответственно.

Деформации, согласно формуле (1), считают малыми, однако перемещения узлов при этом могут быть сколь угодно большими. Равномерно распределенную вдоль участка стержня массу заменяем сосредоточенными массами та, тъ в узлах КЭ (см. рис. 2) [12]:

та = тъ =

рЛ/0

(2)

Векторные уравнения движения составляем на основании второго закона Ньютона, согласно схеме нагрузок, действующих на узел КЭ-модели (рис. 3)

й 2г .■

т

Ж2

• = Nм + N ■ + Р.- + т* (3)

где ■ — номер узла КЭ-модели; т- — узловая масса, состоящая из сумм узловых масс соседних КЭ;

Рис. 3. Усилия, действующие на узел КЭ-модели

г у — радиус-вектор узлау, гу = (х, у, г)т ; Nу, N7_1 — силы, действующие со стороны КЭ с номерами у и (г,-) у - 1, сходящихся в у-й узловой точке; Ру — вектор внешних сил, действующих на у-ю узловую точку; ту g — сила тяжести, действующая на у-ю узловую массу.

Уравнения вида (3) были записаны для всех подвижных узлов КЭ-модели и дополнены начальными условиями. Таким образом, расчет стержня оказался сведенным к решению задачи Коши для системы обыкновенных дифференциальных уравнений (4) (см. ниже) относительно векторов г у (V):

т

й 2гу у й2"

= N у_1 + N у + Ру + туg,

йг у / \

ИТ ('0 ) =у

у 0'

(4)

гу (0 ) = у = 2, ..., п,

у о,

где (0 — начальный момент времени; п — количество КЭ в модели.

Векторы г у (V), определяющие положение узловых точек в пространстве, и будут решением системы (4) с учетом граничных условий г10) = г1 (?о), гп+1 0) = гп+1 (^0).

Определение текущих усилий в стержнях. В уравнения (4) входят текущие усилия в стержнях, которые, согласно (1), связаны с изменением длины КЭ и, следовательно, с узловыми перемещениями. Поскольку в данной постановке текущие усилия в стержнях — векторные величины, а в качестве неизвестных приняты радиусы-векторы узловых точек, то целесообразно определить направляющий орт для каждого элемента и представить изменение длины стержня через соответствующие радиусы-векторы.

Изменение длины элемента с номером у зависит от положения узлов в момент времени I (рис. 4) и исходной длины стержня:

М = I _ ¡0 = |Дг| _ ¡0; (5)

Лг = гу + 1 _ гу . (6)

Тогда усилие ву-м стержне определяем в соответствии с (1) и (5):

N = Ые = —ЕАе, (7)

¡0

е =

Лг

|Лг|

(8)

где е — направляющий орт у-го элемента.

Рис. 4. Положение узлов КЭ и усилия в стержнях в деформированном состоянии

Учет ветровой нагрузки. В задаче о колебаниях и динамической устойчивости электрического кабеля в потоке воздуха (рис. 5) математическая модель должна содержать адекватное описание ветровых нагрузок.

Рис. 5. Электрический кабель в потоке воздуха

Для получения аналитических выражений аэродинамических сил, действующих на стержень, движущийся в потоке, как правило, применяют экспериментально-теоретический метод, при этом полагают, что аэродинамические силы зависят от относительной скорости [2, 13]:

!' (9)

где V — скорость потока; и нии стержня.

вектор перемещений точки осевой ли-

В [2] получены следующие аналитические выражения для распределенных аэродинамических нагрузок:

М = 2 С1Р0й (VI от )2,

1 (10) Ы = 2СпР0й (п от )2,

где — векторы касательной и нормальной распределенной

нагрузки; с1, сп — аэродинамические коэффициенты, определяемые экспериментально; р0 — плотность обтекающей среды; й — характерный размер сечения; v1 от, vn от — касательная и нормальная составляющие относительной скорости.

В общем случае могут присутствовать и другие распределенные нагрузки, например сила Кармана [14], волновое сопротивление [15], которые в данной работе не приведены, так как рассматривается безотрывное обтекание однородным потоком.

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

V■ + V■+1

V- от = V--, (11)

2

где V., V)■ + 1 — скорости узловых точек.

Рис. 6. Замена распределенной аэродинамической нагрузки сосредоточенными узловыми силами

Данное предположение тем точнее, чем меньше длина КЭ. Узловую нагрузку, приложенную в --м и (- + 1)-м узлах, от распределенных аэродинамических нагрузок, действующих на --й элемент, определяем следующим образом:

Fj = Fj+i = ^ = (Ыe + |q„|n)) , (12)

где e, n — орты, указывающие направления q1 и qn. Орт e для касательной нагрузки определен соотношением (8). Орт n для нормальной нагрузки обусловлен исключением из полной скорости касательной составляющей:

v /йот v j от _ v j1 от

n =^-Г = ^-—|. (13)

v jn от v j от — v j 1 от

Узловые аэродинамические нагрузки от соседних элементов суммируем в узлах, соединяющих элементы.

По описанному методу был решен ряд задач, в том числе несколько тестовых. Алгоритм разработан при использовании компьютерного пакета Wolfram Mathematica 9.0.1 [16].

Расчет формы равновесия троса. Очевидно, что метод моделирования движения должен правильно определять и положение равновесия. Тестовая задача о равновесии растяжимого троса рассмотрена в [11].

Концы троса закреплены на опорах на одном уровне. Пролет составляет 1000 м, длина недеформированного троса равна пролету, что характерно для реальных тросовых систем, EA = 1,310 • 107 Н — жесткость сечения на растяжение, pA = 1 кг/м — распределенная масса. Материал троса линейно-упругий.

Количество элементов принято равным 200. Тогда общее число узлов составит 201, число подвижных узлов — 199. Трос отпускали из горизонтального положения без начальной скорости. Для получения равновесной формы троса методом установления в систему уравнений были включены силы сопротивления, пропорциональные массам и скоростям отдельных узловых точек. Внешнее демпфирование, возникающее при отпускании троса, приводило к затухающим колебаниям, после которых трос достигал положения равновесия.

Для сопоставления в таблице приведены результаты вычисления по предложенному алгоритму и решения, полученного в [11].

Результаты решения тестовой задачи двумя способами

Полученные результаты Способ решения Погрешность, %

[11, с. 91] Предложенный алгоритм

Растягивающее усилие в тросе на опоре, Н 37 670 37 657 0,03

Растягивающее усилие в середине пролета, Н 37 347 37 339 0,021

Стрела в середине пролета, м 32,7800 32,7814 0,004

Из таблицы следует, что вычисления, полученные с помощью алгоритма, практически полностью совпадают с результатами, приведенными в [11].

Определение критической скорости потока. Исследовали движение гибкого стержня (кабеля), нагруженного собственным весом и обдуваемого потоком воздуха (см. рис. 5). При некоторой скорости потока могут возникнуть интенсивные вибрации, приводящие к обрыву кабеля.

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

Рассмотрим кабель круглого поперечного сечения диаметром й = 0,1 м, для которого с1 и сп — константы (стоит отметить, что в действительности это справедливо лишь в некотором диапазоне чисел Рейнольдса). В соответствии с [17] принимаем сп = 1,2, с1 = 0,2.

Материал кабеля вязкоупругий, уравнение состояния КЭ:

N =

ЕА

( , йЫ ^

М + , (14)

)

й = е(( -), (15)

где Е = 63 • 108 Па; п = 10-4 с — коэффициент внутреннего трения по гипотезе Фойгта; р — плотность материала кабеля, р = 800 кг/м3. Поток ветра принимали произвольно направленным:

V = с (I + Vyj + к), (15)

где С — безразмерный численный множитель, моделирующий увеличение скорости потока при сохранении его направления; Vх, Vy, vz — проекции вектора V на координатные оси.

Расчет проводили при следующих значениях проекций скорости: Vх = 40 м/с, Vy = 50 м/с, ^ = 50 м/с. Плотность обтекающей среды р0 = = 1,3 кг/м3.

Концы кабеля неподвижно закреплены на опорах на одном уровне. Длина пролета 10 м, длина недеформированного кабеля 24,2 м. Количество элементов принято п = 40, тогда общее число узлов равно 41, число подвижных узлов — 39.

Рассматривали временной интервал 0...60 с. Моделирование проводили в два этапа. На первом этапе методом установления определяли статическое положение равновесия кабеля. На втором рассматривали изменение формы кабеля в случае приложения ветровой нагрузки.

Разделение на этапы условно, поскольку численное интегрирование проходило непрерывно, а «включение» ветровой нагрузки при г = 3,0 с обеспечивала функция Хевисайда.

Положения кабеля в плоскости ХУ и Х2 в моменты времени г = = 2,0; 3,2; 4,4; 58,8; 59,4; 60,0 с показаны на рис. 7.

Рис. 7. Положение кабеля под действием ветровой нагрузки в плоскости ХУ (вверху) и Х2 (внизу) при С = 0,41:

1 — г = 60,0 с; 2 — г = 59,4 с; 3 — г = 4,4 с; 4 — г = 3,2 с; 5 — г = 2,0 с; 6 — г = 58,8 с

Далее модуль скорости ветрового потока с сохранением направления увеличили и определили значение С = С*, при котором возникали интенсивные вибрации.

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

Положения кабеля при С* = 0,42 в моменты времени г = 58,8; 59,4; 60,0 с, для которых при С = 0,41 конфигурация кабеля была одной и той же, показаны на рис. 8.

Графики наглядно показывают, что при значении С = С*, которое соответствует скорости потока IV | = 34,1 м/с, в системе возникают

К м

ч -1 2 \ X,

з Л 3

-А 1

2, м

8 2 /

6

/ 4 /

2 /

. 1 = 59,4 ■ 1 = 58.8

л-1 = 60.0

. I - 59,4 • 1 = 58,8

.V1 = 60,0

-4 -2 0 2 4 6 X, и

Рис. 8. Положение кабеля под действием ветровой нагрузки в плоскости XX (вверху) и Х2 (внизу) при С = С*:

1 — г = 59,4 с; 2 — г = 58,8 с; 3 — г = 60,0 с

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

Заключение. В рамках работы был реализован алгоритм, позволяющий моделировать динамику гибких пространственных стержневых конструкций. Суть подхода — в составлении уравнений движения для отдельных узловых точек системы после замены континуальной модели на дискретную. К основным преимуществам такого подхода можно отнести: использование обыкновенных дифференциальных уравнений вместо дифференциальных уравнений в частных производных; простоту учета геометрической и физической нелинейности; простоту анализа влияния различных параметров на поведение системы.

ЛИТЕРАТУРА

[1] Светлицкий В.А. Механика абсолютно гибких стержней. Москва, Изд-во МАИ, 2001, 432 с.

[2] Светлицкий В.А. Механика стержней. В 2 ч. Ч. 2: Динамика. Москва, Высш. шк., 1987, 304 с.

[3] Сорокин Ф.Д. Прямое тензорное представление уравнений больших перемещений гибкого стержня с использованием вектора конечного поворота. Изв. РАН, МТТ, 1994, № 1, с. 164-168.

[4] Шклярчук Ф.Н., Данилин А.Н. Нелинейные колебания и галопирование провода с обледенением. Изв. ТулГу, Сер. Технические науки, 2013, вып. 11, с. 188-197.

[5] Wang X., Lou W.J. Numerical Approach to Galloping of Conductor. Proc. of the 7th Asia-Pacific. Conference on Wind Engineering. Taipei, Taiwan, 2009, 8 p.

[6] Иванова О.А. О выборе базиса для моделирования движения провода ЛЭП методом Галеркина. Наука и образование, 2013, № 9. URL: http:// technomag.bmstu.ru/doc/602290.html

[7] Ванько В.И., Марчевский И.К. Пляска проводов ЛЭП — неустойчивость по Ляпунову. Энергетика. Известия высших учебных заведений и энергетических объединений СНГ, 2014, № 6, с. 14-23.

[8] Shehata A.Y., El Damatty A.A. Failure analysis of a transmission tower during a microburst. Wind and Structures, 2008, no. 11 (3), pp. 193-208.

[9] Shehata A.Y., El Damatty AA., Savory E. Finite element modeling of transmission line under downburst wind loading. Finite Elem. Anal. Des., 2005, no. 42, pp. 71-89.

[10] Hamada A., El Damatty A.A. Behaviour of guyed transmission line structures under tornado wind loading. Computers and Structures, 2011, no. 89 (11-12), pp. 986-1003.

[11] Зылев В.Б. Вычислительные методы в нелинейной механике конструкций. Москва, НИЦ «Инженер», 1999, 145 с.

[12] Wu S.R. Lumped mass matrix in explicit finite element method for transient dynamics of elasticity. Comput. Methods Appl. Mech. Eng., 2006, no. 195, pp. 5983-5994.

[13] Lilien J.L., Snegovski D. Hurricane Simulation on Power Transmission Line. Proc. Vth Cable Dynamics Symp, Santa Margherita, 2003, pp. 313-318.

[14] Marchevskii I.K., Ivanova O.A. Numerical simulation of wind resonance of a circular profile by means of the vortex element method. J. of Machinery Manufacture and Reliability, 2009, vol. 38, no. 5, pp. 420-424.

[15] Владимиров И.Ю., Корчагин Н.Н., Савин А.С. Гидродинамические реакции при циркуляционном обтекании трубопровода придонным морским течением. Математическое моделирование и численные методы, 2015, № 3, с. 41-57.

[16] Дьяконов В.П. Mathematica 5.1/5.2/6. Программирование и математические вычисления. Москва, ДМК-Пресс, 2008, 574 с.

[17] Соколов А.И. Нелинейные колебания абсолютно гибкого провода в потоке воздуха. Наука и образование, 2008, № 4. URL: http://technomag.edu.ru/ doc/87224.html

Статья поступила в редакцию 25.02.2016

Ссылку на эту статью просим оформлять следующим образом:

Сорокин Ф.Д., Низаметдинов Ф.Р. Численное моделирование движения абсолютно гибкого стержня в потоке воздуха. Математическое моделирование и численные методы, 2016, № 1 (9), с. 3-16.

Сорокин Федор Дмитриевич — д-р техн. наук, профессор кафедры «Прикладная механика» МГТУ им. Н.Э. Баумана. Автор более 60 публикаций в отечественных и международных научных журналах. e-mail: [email protected]

Низаметдинов Фярит Ринатович — студент кафедры «Прикладная механика» МГТУ им. Н.Э. Баумана. e-mail: [email protected]

Numerical simulation of absolutely flexible bAR motion

in the air flow

©F.D. Sorokin, F.R. Nizametdinov Bauman Moscow State Technical University, Moscow, 105005, Russia

The article offers the calculation algorithm of deflected mode of an absolutely flexible bar interacting with the external air flow. The algorithm is based on the replacement of the continual mechanical system by the discrete set of rectilinear finite elements and concentrated masses. The authors show differential equations of mass motion with allowance for an aerodynamic load and dissipative forces and integrate them by numerical method. That made it possible to find both the equilibrium position of the flexible bar in the flow, and the critical flow velocity which causes violent bar vibrations in case of its excess.

Keywords: absolutely flexible bar, aerodynamic load, finite element, numerical integration, dynamic instability, autooscillation.

REFERENCES

[1] Svetlitskyy V.A. Mekhanika absolyutno gibkikh sterzhney [Mechanics of absolutely flexible bars]. Moscow, Moscow Aviation Institute Publ., 2001, 432 p.

[2] Svetlitskyy V.A. Mekhanika sterzhney. Dinamika [Mechanics of bars. Dynamics]. Moscow, Higher school Publ., 1987, vol. 2, 304 p.

[3] Sorokin F.D. Izvestiya RAN — Proceedings of RAS. MTT, 1994, no. 1, pp. 164-168.

[4] Shklyarchuk F.N., DanilinA. N. Izvestiya Tulskogo Gosudarstvennogo universi-teta. Tekhnicheskiye nauki — Proceedings of Tula State University.Technical Science, 2013, vol. 11, pp. 188-197.

[5] Wang X., Lou W.J. Numerical Approach to Galloping of Conductor. Proc. of the 7-th Asia-Pacific. Conference on Wind Engineering. Taipei, Taiwan, 2009, 8 p.

[6] Ivanova O. A. Nauka i obrazovaniye. Elektronnyy Zhurnal — Science and Education. Electronic journal, 2013, no. 9. Available at: http:// tech-nomag.bmstu.ru/doc/602290.html

[7] Vanko V.I., Marchevskyy I. K., Izvestiya visshikh uchebnikh zavedenyy i ener-geticheskikh obyedinenyy SNG — Proceedings of institution of higher education and power associations of the CIS, 2014, no. 6, pp. 14-23.

[8] Shehata, A. Y., El Damatty, A. A. Wind and Structures, 2008, no. 11(3), pp. 193-208.

[9] Shehata, A.Y., El Damatty, A.A., Savory, E . Finite Elem. Anal. Des., 2005, no. 42, pp. 71-89.

[10] Hamada, A., and El Damatty, A. A Computers and Structures, 2011, no. 89 (11-12), pp. 986-1003.

[11] Zylyov V.B. Vychislitelniye metody v nelineynoy mekhanike konstruktsyy [Computational approach in nonlinear mechanics of constructions]. Moscow, NITS Inzhener, 1999, 145 p.

[12] Wu S.R. Comput. Methods Appl. Mech. Engrg, 2006, no. 195, pp. 5983-5994.

[13] Lilien J.L., Snegovski D. Proc. Vth Cable Dynamics Symp., Santa Margherita, 2003, pp. 313-318.

[14] Marchevskii I.K., Ivanova O.A. Journal of Machinery Manufacture and Reliability, 2009, vol. 38, no. 5, pp. 420-424.

[15] Vladimirov I.Y., Korchagin N.N., Savin, A.S. Matematicheskoye modeliro-vaniye i chislennyye metody — Mathematical Modeling and Computational Methods, 2015, no. 3, pp. 41-57.

[16] DyakonovV.P. Programmirovaniye i matematicheskiy evychisleniya [Programming and mathematical calculations]. Moscow, DMK-Press Publ., 2008, 574 p.

[17] Sokolov A. I. Nauka i obrazovaniye. Elektronnyy Zhurnal — Science and Education. Electronic journal, 2008, no. 4. Available at: http://technomag.edu.ru/ doc/87224.html

Sorokin F.D. — Dr. Sci. (Eng.), Professor of the Department of Applied Mechanics, Bauman Moscow State Technical University. Author of more than 60 publications in the national and international magazines. e-mail: [email protected]

Nizametdinov F.R. — a student of the Department of Applied Mechanics, Bauman Moscow State Technical University. e-mail: [email protected]