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

УДК 004.94+629.78

Создание

математической модели общепланетарного транспортного средства на основе метода конечных элементов

А.Э. Юницкий12

доктор философии транспорта

А.Р. Покульницкий2

1 ООО «Астроинженерные технологии»,

г. Минск, Беларусь

2 ЗАО «Струнные технологии», г. Минск, Беларусь

ГГ

Описаны подходы к разработке программного комплекса, предназначенного для математического моделирования поведения общепланетарного транспортного средства (ОТС) с помощью универсальных расчётных элементов. В основу положен метод конечных элементов. Рассмотрен процесс создания балочного конечного элемента, подходящего для характеристики поведения ОТС. Учтены нелинейности летательного аппарата, обусловленные его сложной геометрией, нерегулярной структурой и особенностями нагружения. Приведены библиотеки, которые применялись для реализации программного комплекса, аргументировано их использование исходя из модульности продукта и отдельных функций каждого модуля. Осуществлена верификация конечного элемента на основе базовых задач механики и представлено последующее сравнение с альтернативными вариантами решения. Выполнена апробация элемента на задачах Струнного транспорта Юницкого (ЮСТ) ввиду схожести их нелинейности, а также достаточного количества экспериментальных данных для оценки точности, устойчивости и сходимости разработанного конечного элемента.

Ключевые слова:

балочный конечный элемент, математическая модель общепланетарного транспортного средства (ОТС), метод конечных элементов, метод Ньютона - Рафсона, элемент Эйлера - Бернулли, язык программирования Python.

UDC 004.94+629.78

Creation

of the Mathematical Model of the General Planetary Vehicle on the Basis of the Method of Final Elements

A. Unitsky12

Ph.D. in Information Technologies (Transport)

A. Pakulnitski2

1 Astroengineering Technologies LLC, Minsk, Belarus

2 Unitsky String Technologies Inc., Minsk, Belarus

PP

The authors describe approaches to the development of a software complex designed for mathematical modelling of the General Planetary Vehicle (GPV) behavior with the help of common calculation elements. Such approaches are based on the finite element method. There is a review of the process of creation of a beam-type finite element that can be used for the GPV behavior description. The aircraft nonlinearities arising from its complex geometry, irregular structure and loading features are considered. The article presents the libraries which were used in the course of the software complex designing, as well as the justification of their use according to the modularity of the product and individual functions of each module. Verification of the finite element based on the core issues of mechanics and further comparison with the alternative solutions have been performed. The element has been tested on tasks of the Unitsky String Transport (uST) due to the similarity of their nonlinearity as well as sufficient experimental data to assess the accuracy, stability and convergence of the developed finite element.

Keywords:

beam-type finite element, Euler-Bernouili element, finite element method, mathematical model of the General Planetary Vehicle (GPV), Newton-Raphson method, Python programming language.

Введение

Создание математической модели общепланетарного транспортного средства (ОТС) [1-3] - основополагающий этап при выполнении реальных инженерных расчётов в данной области. С учётом сложности конструкции, большой длины (40 076 км), огромной массы (до 40 млн тонн), высокой скорости движения ленточных маховиков (до 12 км/с) и различных режимов работы летательного аппарата, а также ввиду невозможности проведения целого ряда натурных испытаний математическое моделирование ускоряет работу, экономит ресурсы и обеспечивает вариативность решений такого рода прикладных задач. В частности, оно позволяет описывать режимы работы и процессы, происходящие с ОТС на стадии взлёта, функционирования на экваториальной орбите и во время посадки, а также прогнозировать нештатные ситуации, которые могут возникнуть во время эксплуатации, и планировать возможные пути их недопущения.

В [1, 3] изложен аналитический метод получения кинематических характеристик движения ОТС, рассмотрены вопросы, связанные с затратой энергии при его выведении на орбиту и эксплуатации. Настоящее исследование описывает процесс создания программного комплекса, предназначенного для математического моделирования поведения ОТС посредством применения универсальных расчётных элементов балочного типа.

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

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

создавать необходимые модели и расчётные схемы, оптимизировать расчётное время, мощностные затраты, а также решать задачи масштабирования и преемственности.

Основные подходы

к моделированию программного комплекса

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

• подобрать жёсткостно-массовые и жёсткостно-демп-фирующие характеристики элементов конструкции ОТС;

• определить собственные частоты колебаний ОТС;

• провести гармонический анализ и анализ, связанный со случайной вибрацией;

• выполнить анализ напряжённо-деформированного состояния.

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

• неравномерная тяга линейных электродвигателей по длине ОТС;

• неравномерный КПД линейных электродвигателей по длине ОТС и его изменение во времени (разгон и торможение ленточных маховиков, в том числе в зависимости от скорости их движения в диапазоне 0-12 км/с);

• неравномерность массы по длине ОТС (прежде всего из-за неравномерного размещения пассажиров и грузов);

• начальные геометрические отклонения от идеальной окружности (повторение формы геоида на океанических участках и крупного рельефа местности на суше и при пересечении гор);

• неравномерность гравитационного поля планеты (рисунок 1);

• нелинейности сил магнитов в магнитной подвеске линейных роторов (маховиков);

• иное.

Аномалия гравитационного поля Земли, мГал

Рисунок! - Карта аномалий гравитационного поля Земли

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

[\К\ • \К\• [5,] • Ш • М){») = - {/ },(1)

где [Ке] - матрица линейной жёсткости;

[Кд] - матрица геометрической жёсткости; [5,1 [5Д [5,] - матрицы, позволяющие учесть в элементе эффекты деформации более высокого порядка малости. В свою очередь, матрицы получены в следующем

виде:

- ЕА I 0 0 ЕА I 0 0

0 12 Е1 Е 60 Е 0 120 Е 60 Е

Кг- 0 ЕА I 60 Е 0 40 0 0 ЕА I 60 " Е 0 20 0

0 120 Е 60 " Е 0 120 Е 60 " Е

0 60 Е 20 1 0 60 " Е 40 I _

В основу разрабатываемого программного комплекса положен балочный элемент Эйлера - Бернулли [4], внешний вид которого представлен на рисунке 2. Данный элемент описывает линейное поведение, однако с учётом приведённых выше нелинейностей в математическую модель добавлены эффекты деформации более высокого порядка малости и геометрическая жёсткость. В зависимости от применяемых подходов к описанию уравнений и их решению создаваемый элемент универсален и позволяет моделировать задачи как статики, так и динамики.

Рисунок 2 - Балочный конечный элемент: а - перемещение и углы поворотов узлов; б - силовые факторы в узлах

' хЬ

~г

0 ^

1

!~хь

м*

0 М ' 2Э Я*

I I

6^-6 12 Я* 0

51 А1} 10 1 Ж2 51

Я* Ш}:Ь Щь '•// м /3 Я*

10 АС 15 ' М_ I 10

0 М ' ¿3 Я*

I

12 1С,, 0

51 10 Ж2 51

Я* ^хь 2//^, Я*

10 ' Ж2 30 1 10

к.

АС АС

А15

А1}

[

я* щь

тг АС

Щь

сэ А1

М*

I

^ б/я.

10 АС

ЩЬ ^хь

15

А1 J

М

гап ^~хЬп ^~уЬп

(4)

' хЫ

~г

1~хЫ

м„

' хЫ

~г

м„

1~хЫ

0 12/^ Я« 6/^ 0 12

51 /II5 10 АС 51 + /II5

Я'Ы Щы Щы Я'Ы 6/^-ы

10 + Л/.2 30 ' Ж 10 + ж2

к.

0 Щы 0 Щы

51 + АС 10 1 Ж2 51 АС 10 1 АС

Я'Ы 6/^-6/ У-Щы Щы Щы

10 Ж2 15 ' А1 1 10 АС 30 1 А1

к

т

1~хЫ Щы

10 АС

Щь,

15 ' Ж

' хЫ

~г

' хЫ

~г

м„

^ О

м7а

' хЫ

~г

0 6 СхЫ Я« 0 Щы 1~хЫ

51 + АС 10 1 АС 51 АС 10 1 АС

61СхЬ! У-Щы ^ы Я'Ы Щы

1 10 АС 15 ' Ж 1 10 АС 30 ' Ж

и*

6 Схы 12 ^хЫ ^хы 0 ^хы 12 ^хЫ

51 Ж5 10 АС 51 ' Ж5 10

Я« 6/^-6/ 1~Схы щы Я'Ы

10 1 АС 30 ' Ж 10 1 И12 15 1 ж

где Е- модуль упругости;

А - площадь поперечного сечения; ¿-длина элемента; /- момент инерции;

Ехап' Еузп' Ехьп' Еуьр, /"в/, Е г Еу^, М , М , М , М -

компоненты сил и моментов, входящие в состав матриц (2)-(6) (описаны в (7), (8));

индексы а и Ь - приложение нагрузки в узле а или Ь (рисунок 2); х,у, I- оси координат; п и /- нелинейные и линейные компоненты сил.

хз1 СН £

Еы= ~Еха1,

Еуы Еуд1,

Е!

(7)

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

Метод обладает квадратичной сходимостью и для решателя программного комплекса может быть представлен в виде:

Г = ^ - Е!

гег е\1. /г/Г'

р < критерий сходимости,

(9) (Ю) (11)

где Ди" - приращение перемещения на итерации; приложенная нагрузка; /^г- рассчитанное внутреннее усилие; Кп - жёсткость системы;

разница двух усилий (критерий по силе). Динамическое поведение элемента описывается основным уравнением динамики, приведённым в матричном виде:

Г = --ЕА

1 хап 2

и, - и,

ЕА_

2 Е

5а 15

1 ' 1 5 % + 15 ч, ■

+ 1 - ^К + - - ^ШЕ - -¿Ий/

12 1 + 5

15

15

. Ъ;Е. а, . . ш-• Ш-

г

уз/7

м

1 Ьп ¡2

и, - и„

I

и„ - иР

(12и + бе/-12и6 + беД

(8)

^~хЬп ^'хап1 ^~уЬп ~ ^~уап1

м

1 ':Ьп ['2

и, - и,

I

(6|/ + 26/-6|/6 + 4еД

где и3, иь- перемещение узлов вдоль оси Х\ vb - перемещение узлов вдоль оси К; 63,6Ь - углы поворота узлов (рисунок 2). Для решения уравнения (1) в программный комплекс добавлен метод Ньютона - Рафсона в полной постановке [5] -

[ММ + ЮМ + Я» =

(12)

где [М] - матрица масс балочного элемента; [С] - матрица демпфирования. Для решения задач динамики применялся метод Ньюмарка, который может быть представлен следующим образом:

й^ = йп + Щ, йу = [1-у)й„ + уй„^ й„х1 = йп + ['\-у}Шп + уЩ1,],

йр = [1-2р)йп + 2рйп+1, (13) АЁ

V = 0,5, р = 0,25,

где и, и, й - перемещение, скорость, ускорение; ДГ- временной шаг.

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

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

• модуль для создания геометрии и описания свойств материала - для редактирования входной геометрии;

• модуль сеточного предпроцессора - для представления геометрической модели в виде математической;

• модуль приложения нагрузок и граничных условий -для задания условий проведения расчёта;

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

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

F= 5000 Н

Рисунок 3 - Расчётная схема задачи 1. Консольная балка (статическая постановка)

F= 5000 Н

Рисунок 4 - Вертикальный прогиб консольной балки (статическая постановка)

F= 310 Н

Верификация полученных результатов

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

1. Расчёт длинной консольной балки (¿ = 2,5 м) квадратного сечения (50 * 50 мм), нагруженной на конце сосредоточенной силой [F= 5000 Н) (рисунок 3). Материал - сталь, модуль упругости £= 200 000 МПа, коэффициент Пуассона v = 0,3. Результаты моделирования представлены на рисунке 4.

Следует отметить, что задача о консольной балке выбрана для верификации конечного элемента по причине того, что соотношение высоты к длине балки составляет 1/50 и имеет нелинейное поведение (большие деформации). На рисунке 4 показана деформация балки под действием заданной силы. Максимальный прогиб на конце балки составил 247,53 мм.

2. Расчёт конструкции в виде двух стержней, жёстко защемлённых по концам. В средней точке стержни соединены с подъёмом 9,8 мм. Геометрические размеры представлены на рисунке 5. Материал стержней - алюминий, модуль упругости Е= 70 922 МПа, коэффициент Пуассона v = 0,2, нагрузка F= = 310 Н. Результаты моделирования представлены на рисунке 6.

STTi

Рисунок 5 - Расчётная схема задачи 2, геометрические размеры образца. Классическая задача с распором (статическая постановка)

F= 310 Н

Z = 19,02 мм

Рисунок 6 - Вертикальный прогиб конструкции из двух стержней (статическая постановка)

Такая задача важна тем, что в процессе деформирования происходит нелинейное изменение жёсткости конструкции и изменение усилий в стержнях (распор). На рисунке 6 показана деформация балок под действием заданной силы. Максимальный прогиб в середине пролёта составил 19,02 мм.

3. Для проверки корректности динамической составляющей решения применялась задача о колебаниях консольной балки (динамическая постановка) квадратного сечения (50 * 50 мм), нагруженной на конце сосредоточенной постоянной силой [F= 5000 Н) (рисунок 7). Материал -сталь, модуль упругости Е = 200 000 МПа, конструкционный коэффициент затухания = 0,01. На рисунках 8 и 9 представлены результаты моделирования, а также динамическое поведение консольной балки, нагруженной постоянной силой. Рисунок 9 - графическое изображение колебаний свободного конца балки. Как видно, максимальный прогиб составляет 376,05 мм. После затухания его величина совпадает с полученной в статической задаче 1.

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

0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9

I

F= 5000 Н

50 мм

J

¿ = 2,5 м

<- -►

-50 -

? -100 -

21

ai" s -150 -

CD =Г -200 -

ш

CD о_ -250 -

1= -300 -

-350 -

-400 -

Максимальный прогиб: 1= 376,05 мм Время, с

Рисунок 9 - Динамический отклик консольной балки, нагруженной постоянной силой

Таблица - Результаты моделирования

Рисунок7 - Расчётная схема задачи 3. Консольная балка (динамическая постановка)

Показатель Номер задачи Разрабатываемый программный продукт ANSYS Аналитическое решение

Вертикальное перемещение, мм 1 247,53 247,84 247,67

2 19,02 19,36 19,27

Динамический прогиб, максимальное значение, мм 3 376,05 376,75 376,32

Динамический прогиб, значение после затухания, мм 3 247,53 247,84 247,67

F= 5000 Н

Z= 247,53 мм

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

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

Апробация полученных результатов

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

Внешний вид поперечного сечения путевой структуры представлен на рисунке 10. Её геометрические и физические характеристики: вес нижнего пояса - 84 кг/пог. м; вес верхнего пояса - 15,5 кг/пог. м; длина - 720 м; натяжение монтажное, нижний пояс - 5730 кН; натяжение монтажное, верхний пояс - 1860 кН; общее натяжение - 7240 кН; шаг подвесов - 3 м.

На рисунке 11 показана упрощённая модель транспорта. Его механические характеристики: общая жёсткость первой ступени подрессоривания - 220 Н/мм; общая жёсткость второй ступени подрессоривания - 520 Н/мм; демпфирование первой ступени подрессоривания - 9 Н-с/мм; скорость проезда - 60-100 км/ч; первая собственная частота подвески -1,05 Гц; масса - 8-24 тонн.

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

Рисунок 11 - Упрощённая модель транспортного средства (юнимобиля)

Базовый объём задач в рассматриваемом случае связан с подбором траектории движения юнимобиля, при которой в пассажирском модуле возникают вертикальные ускорения, комфортные для проезда людей. Результаты моделирования представлены на рисунках 12 и 13.

В качестве практического применения программного комплекса для задач, связанных с моделированием поведения ОТС, проведены кинематические расчёты перемещения, скорости и ускорения сегментарного участка данного летательного аппарата (рисунок 14).

Полученные с использованием метода конечных элементов характеристики (рисунок 15) полностью соответствуют вычисленным ранее А.Э. Юницким аналитическим методом [1,2].

гГЬ, гНп

Рисунок 10 - Сечение путевой структуры

Промежуточная опора 1 Промежуточная опора 2

5 -4 -3 -

ъ 2 -

0 -

-2 -I

-3 -4

Ускорение при проезде опоры: -0,7 м/с;

Ускорение при проезде опоры: -0,3 м/сг

60

65

70

75

Время, с

85

90

Рисунок 12 - Динамический отклик в пассажирском модуле (масса юнимобиля 8000 кг, скорость движения 100 км/ч; вертикальная составляющая)

Рисунок 13 - Полученная оптимальная траектория движения транспорта

2000 3000 Время, с

Длина, м

— Верхний пояс Нижний пояс

Анкерная опора 2

0 100 200 300 400 500 600 700

800

186 -Анкерная опора 1

0 1000 2000 3000 4000 5000

Время, с

Рисунок 14 - Расчётная схема для определения кинематических характеристик (динамическая постановка задачи): V- взлётные скорости компонентов ОТС; /?- радиус Земли; в- сила притяжения Земли; О - сила сопротивления атмосферы;

Г- сила упругости ОТС; ы - угловая скорость Земли

Заключение

Разработанный балочный конечный элемент на основе элемента Эйлера - Бернулли при дополнении его функционалом геометрической матрицы жёсткости и матриц учёта деформаций большего порядка малости даёт возможность учитывать геометрическую нелинейность (большие деформации), что позволяет моделировать нелинейности, влияющие на поведение ОТС.

Результаты тестовых задач сравнивались как с полученными аналитическим (точным) методом, так и с вычисленными коммерческим конечно-элементным продуктом ДЫБУЗ.

Рисунок 15 - Полученные кинематические характеристики: а - вертикальное перемещение сегмента ОТС; б - вертикальная скорость сегмента ОТС; в - вертикальное ускорение сегмента ОТС

Расчётные задачи демонстрируют, что созданный конечный элемент даёт хорошую точность и устойчивость решения.

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

В дальнейшем с помощью рассматриваемого программного комплекса будут оптимизированы все требования, предъявляемые к элементам ОТС, для эффективного

и безопасного функционирования этого гигантского геокосмического летательного аппарата: допустимая неравномерность тяги линейных электродвигателей; допустимая неравномерность загрузки ОТС по длине пассажирами и грузами; допустимая неравномерность КПД работы линейных электродвигателей по своей длине и во времени; допустимые горизонтальные и вертикальные отклонения корпуса ОТС (и движущихся внутри с космическими скоростями линейных маховиков) от формы идеальной окружности и др.

Список основных источников

1. Юницкий, A3. Создание математической модели общепланетарного транспортного средства с оптимизацией положения и учётом влияния динамических факторов / А.Э. Юницкий, Р.А. Шаршов, А.А. Абакумов //Безракетная индустриализация ближнего космоса: проблемы, идеи, проекты: материалы III междунар. науч.-техн. конф, Марьина Горка, 12 сент. 2020 г. / ООО «Астроинженерные технологии», ЗАО «Струнные технологии»; под общ. ред. A3. Юницкого. -Минск:СтройМедиаПроец2021.-С. 182-193.

2. Юницкий, A3. Создание математической модели общепланетарного транспортного средства: разгон маховиков, прохождение атмосферы, выход на орбиту/A3. Юницкий, Р. А. Шаршов, АЛ Абакумов // Безракетная индустриализация космоса: проблемы, идеи, проекты: материалы II междунар. науч.-техн. конф., Марьина Горка, 21 июня 2019 г. / ООО «Астроинженерные технологии»; под общ. ред. A3. Юницкого. - Минск: Парадокс, 2019. - С. 77-83.

3. Юницкий, A3. Струнные транспортные системы: на Земле и в Космосе: науч. издание / A3. Юницкий. - Минск: Беларус. навука, 2017. - 342 е.: ил.

4. Kim, N.-H. Introduction to Nonlinear Finite Element Analysis / N.-H. Kim. - New York: Springer, 2018. - 437p.

5. Reddy, J.N. An Introduction to Nonlinear Finite Element Analysis: With Applications to Heat Transfer, Fluid Mechanics and Solid Mechanics /J.N. Reddy. - 2nd ed. - New York: Oxford University Press, 2015. - 768 p.