Исследование нового осесимметричного конечного элемента для моделирования упруго-динамических свойств газа Текст научной статьи по специальности «Физика»

МАШИНОСТРОЕНИЕ И МАШИНОВЕДЕНИЕ ОМСКИЙ НАУЧНЫЙ ВЕСТНИК № 1 (107) 2012

150

И. А. ТРИБЕЛЬСКИЙ В. В. БОХАН

Научно-производственное предприятие «Прогресс», г. Омск

ИССЛЕДОВАНИЕ НОВОГО ОСЕСИММЕТРИЧНОГО КОНЕЧНОГО ЭЛЕМЕНТА ДЛЯ МОДЕЛИРОВАНИЯ УПРУГО-ДИНАМИЧЕСКИХ СВОЙСТВ ГАЗА________________________________________________________

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

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

В [1] описан новый конечный элемент (КЭ), моделирующий газ — его упругие и инерционные свойства как в динамических задачах, так и в статических.

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

В настоящее время разработан ряд коммерческих программ для расчёта методом конечных элементов (МКЭ), таких как Амуз, На81;гап, Сс^тс^/БоШ-Шогкв, АЪадш. Большинство из них распределение давления при малых колебаниях газа и жидкости описывают волновым уравнением, которое составлено относительно давления. Решением также является давление и, соответственно, КЭ имеет одну степень свободы — давление. Кроме того, в связанных задачах для КЭ на границе твердого тела совместно с волновым дополнительно необходимо решать задачу деформирования твердого тела, описываемую конечными элементами другого типа, что вызывает определённые трудности при стыковке КЭ разных типов.

Новый КЭ позволяет преодолеть указанные ограничения. Он прост, адекватен в описании газа, детально проработан. Для практического применения используется специально разработанная в ФГУП «НПП "Прогресс"» программа. В настоящей работе показано, что новый КЭ можно моделировать в различных программах, реализующих МКЭ, применительно к осесимметричной задаче. Это достигается введением некоторых правил при составлении матрицы жёсткости КЭ. Получено в явном виде выражение для матрицы жёсткости, записанное относительно координат узлов конечного элемента. Таким образом, предложенный в [1] подход, находит более широкое применение.

І і

К гм/ /

/

И' 4

т

Рис. 1. Треугольный осесимметричный КЭ цт

Рассмотрим вывод матрицы жёсткости треугольного осесимметричного КЭ двумя способами: способом, который предложил И. А. Трибельский в [1], и способом, описанным в работах [2 — 5]. Последний способ широко используется в программах для расчёта прочности, устойчивости, отклика на воздействие различных конструкций.

На рис. 1 изображён треугольный осесимметричный КЭ с узлами і, ], т в системе координат г, z.

Рассмотрим вывод выражения матрицы жёсткости КЭ способом, описанным в работе [1].

Запишем объём V тетраэдрального КЭ через координаты узлов:

1 г- гі

і 1 г. z,

2 1 ]

гт 2т

= 2%г

1 (г

2

+ г,2, - г,2, )

(1)

т-а

т-а

УДК 62-784.3

%

Здесь 5 — площадь треугольника у'т, гтШ=(г+ +г]+гт)/3 — средний радиус КЭ, г, г, гш, 2,, 2, zш — координаты узлов.

Изменение объёма КЭ запишем через изменения координат его узлов:

Атг ЭУ„ ЭУ„ ЭУ„ ЭУ ,

А У = --- ОГ. +--02; +------0Г; +----02 ; +

Эг; Э2; Эгц 1 Э2ц 1

ЭУ , ЭУ ,

+ ^°Гш +^&ш Эг Э2„

- (Ро + Ap)d(AV) = -= (й{8})т ^ ]

= №})

V

Приравнивая (10) и (11), получим:

(11)

или в матричном виде

ДV=[DV]{8}.

(2)

Здесь [ БУ] — вектор-строка производных объёма КЭ по координатам

{Р} = Е°П ^]Т ^]{8} - Ро [DV]т .

V

(12)

Здесь первый член выражает упругую силу, а второй обусловлен номинальной нагрузкой. Таким образом, из получаем матрицу жёсткости КЭ:

эу _ду

дг, Эг, Эг, дг} Эгт

(3)

{8} — вектор-столбец перемещений по координатам

{8}Г=[8г, 82 Ц 8гш &Ш]. (4)

Подстановка (1) в (3) даёт

[ОУ]=РгтШ[-2т+2 гт-г1 2т-21 -гт+г1 -2+2 —Ь (5)

Всякое изменение объёма КЭ вызывает изменение давление газа. Для политропного процесса можно записать:

Рo(Vo)Л=(Рo + Дp)(Vo + ДУ)Л,

(6)

[К] = Е^П [DV] [DV].

(13)

где р0, V0 — начальные значения давления и объёма газа в КЭ, обусловленные номинальной нагрузкой; Др, ДV — приращения давления и объёма КЭ, вызванные изменением координат узлов; п — показатель политропы.

Разложим в (6) выражение (Vg + Д V)п в ряд:

Выпишем в явном виде зависимость матрицы жёсткости от координат узлов, подставляя (1) и (5) в (13):

[К] = Роп

рг„

25

[К ]ц [К ]12

[К ]21 [К ]22

где для удобства записи обозначено

(- г т + г, )2 (- гт + г, )(гт - г, )

[К]11 =

(гт - г, )(- г т + 2 ) (гт - ^ )2

(гт - гі )(- гт + г, ) (гт - г, )(гт - г, )

(- гт + )(гт - гі )

(гт - г, )(гт - гі )

(гт - г, )2

(14)

, (15)

Vnn

(V0 + AV )п = V; + ^-п

AV +1 Уnnn - ^ 2

AV2 +

1 ^(п - 1)(п - 2)

V

у п

AV3 +.

и удержим в нём первые два члена

Упп

(у +АУ )п = уп + У-уП АУ.

' п

(7)

Подставив (7) в (6), находим (с учётом малости ДрД^:

Ар = - Ео^ AV V

или, с учётом (2),

Ap = - ^ тР}.

* п

(8)

(9)

Матрицу жёсткости [К] получим на основе принципа возможных перемещений. Пусть ¿{8} — некоторое виртуальное перемещение узлов КЭ. Тогда работа внешних узловых сил {Р} будет равна:

И8>т,

(10)

а изменение внутренней энергии КЭ, обусловленное изменением его объёма на величину ¿(ДУ), можно рассчитать с учётом (2) и (9):

[К]12 =

(- гт + 2 )(- гт + гг ) (- гт + г, )(- ^ + г, )

(гт - г, )(- гт + гг )

(гт - г, )(- гт + г■ )

(гт - г, )(- ^ + г, ) (гт - г, )(- 2 + г, )

(- г™ + г,)(г, - Г, ) (гт - г, )(г, - Г, )

(гт - г, )(г, - Г, )

[К]21 =

(- гт + Гі )(- г„ + г; ) (- гш + гг )(ГШ - )

(- г, + г, )(- гт + г, ) (- г, + г, )(гт - г, )

(г, - Ъ)(- гш + г,)

(Г] - г, )(г"т - г, )

(- гт + г, )(гт - гі У

(- 2 + гі )(гт - гі ) (г, - гі )(гт - гі ) .

[К]22 =

(- гт + Т- )2 (- г, + гі )(- гт + г )

(г, - гі )(- гт + гі )

(- гт + г )(- г, + гі ) (- г, + гі )2 (г, - гі )(- 2 + гі )

(- гш + г, )(г, - Г, ) (- г7 + г, )(г, - г, ) (г; - гі )2

(16)

(17)

ОМСКИЙ НАУЧНЫЙ ВЕСТНИК № 1 (107) 2012 МАШИНОСТРОЕНИЕ И МАШИНОВЕДЕНИЕ

МАШИНОСТРОЕНИЕ И МАШИНОВЕДЕНИЕ ОМСКИЙ НАУЧНЫЙ ВЕСТНИК №1 (107) 2012

152

Очевидно, что полученную матрицу жёсткости [К] можно, в принципе, использовать для расчётов с помощью достаточно широкого круга программ, реализующих МКЭ применительно к задачам деформирования твердого тела. Однако далеко не все они позволяют непосредственно ввести значения элементов [ К ]. Зачастую остаётся лишь возможность непосредственного ввода матрицы упругости [Б]. И тем не менее, как показано далее, есть возможность моделировать газ изотропными КЭ аналогично новому КЭ.

Далее рассмотрим вывод выражения матрицы жёсткости КЭ (рис. 1) способом, описанным в работах [2 — 5].

Запишем функции положения:

Несомненно, газ можно принять изотропным материалом, тогда матрица упругости имеет вид:

здесь

1

N. = — (я. + Ьіг + сіг),

і 25 і і 1

N1 = — (а, + Ьг + с ¡г),

і 25 1 1 1

N'т = 25 (ат + Ьтг + стг) ,

а=гг —гг , Ь=—г +г, с=г —г,

і і ш і т і ш 1 і ш 1

а=г г—г г, Ь=г —г, с=—г +г.,

1 Ш і ш і 1 ш і 1 ш і

а =г2—г.г, Ь =—г+г, с =г—г..

ш і 1 і 1 ш 1 і ш 1 і

(19)

(20) (21)

(22)

(23)

(24)

[Б ] =

Е (1 -V)

(1 + п)(1 - 2п)

1 -V 1-V V

1 - V V

V

1 -V 1-V

о о

1 -V 1

о

о

о

о

1 - 2v

2(1 -V)]

, (26)

здесь Е, V — приведённые для газа модуль Юнга и коэффициент Пуассона.

Запишем теперь матрицу жёсткости КЭ:

2п[Б]Т[0][Б]ММ2

[ к]=2р[Б]Т[Б][Б]гшиа5.

(27)

Подставляя в (27) выражения (25) и (26), получим:

, (28)

[к]=ргшіа Е(1 -^

[ к] 11 [к] 12 [к] 21 [к] 22

25 (1 + V)(1 - 2v) где для удобства записи обозначено

Вектор-столбец деформаций для осесимметричной задачи имеет вид:

{е} =

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

е6=0, у =0.

6 ' • г2

Деформации выражают соотношением:

{е} = [Б]{8},

в котором {8} определяется по (4), а матрица [ Б ] определяется выражениями (19) — (24):

[Б ] =

о

N

Эг

о

N

Эг

о

о

о

о ЭN1

Эг

о

[Б] =

25

о

25

о

о

о

ЭN1 Эг о о о

о

гш - гі

25

о

о

N

Эг

о

эNm

Эг

о

о

о

25

о

о

о

■г1 + гі 25

о

25

о

о

о

(25)

[к] и =

(- гш + г1 )2 (- гш + г1 )1^(гт - г1)

(гш - г1 )!^(- гш + г1) (ш - г, )2

(г - г )(- г + г ) (г - г )——(г - г.)

V т і А ш 1 / \ ш і > 1 V ' ш 1 /

(- гш + г1 )(гт - гі )

(гт - г1) ^(гт - гі ) (гт - гі )2

(29)

[к] 12 =

(- гш + г1 )1^(- гш + г ) (- гш + г1 )(- г1 + гі )

(гт - г1 X- гш + г, ) (гт - г1) 1^(- г1 + гі )

(гш - гі ^Т"“ (- гш + г, ) (гш - гі )(- г1 + гі )

1 -V

(- гш + г1 - ^ )

(гт - г1 \г, - г, )

(гш -гі)~г^іг1 -г)

(30)

[к] 21 =

(- гш + г■ )1^(- гш + г1) (- гш + г )(гт - г1)

V 1 -V

(г1 - г- ) 1^(- гш + г1) (г1 - г Ьш - г1)

Г г1 + гі X- гш + г1) (- г1 + гі ):—(гт - г1)

(- гш + г ) 1^(гт - гі )

(- г1 + гі )(гт - гі )

(г1 - г )Т^(гт - ^ )

, (31)

V

V

или

г

о

о

о

г - г.

т1

- гш + г

о

- гш + г1

о

о

г1 - гі

о

о

На рис. 2 показана КЭ-модель трубки диаметром d = 8 см, высотой L = 20 см, заполненной воздухом и с закрытым нижним концом. Противоположный конец трубки — излучатель — перемещается вдоль оси z с амплитудой А =14 см и заданной частотой f. Воздух в трубке находится под давлением p0 = 101325 Па и имеет температуру T = 293,15 К. Параметры газа изменяются по политропному процессу с показателем политропы л =1,4.

Выполнен расчёт амплитуды давления тремя методами: аналитический расчёт; расчёт на основе МКЭ по [1] с использованием специально разработанной в ФГУП «НПП "Прогресс"» программы для нового КЭ; расчёт на основе МКЭ, выполненный в Апзуз с разбиением на изотропные осесимметричные элементы РЬАЫЕ25.

Аналитический расчёт основан на решении волнового уравнения в частных производных [6]:

Э2и 2 Э2и _ , . _

—т = с —т, 0 < г < Ь , і > 0 . ді2 Эг2

(35)

Рис. 2. КЭ-модель закрытой трубки

Здесь и — перемещение частиц среды, с — скорость распространения упругой волны, г — координата вдоль трубки,і —время.

Полагая и=и(г)-віші, при следующих граничных условиях

и(Ь, і)=Асо$(<оі), и(0,і)=0,

[¿1 22 =

(- + Г )2

(- ^ + г )1^(- Гт + Г )

(- г_, + гі )2

(- Гт + Г )(г) - Г )

(- г+г, )т-п(г^- г,)

Ь - Г, )2

, (32)

Сравнивая выражения (14) — (18) с (28) — (32), выявляем условия равенства матриц жёсткости [К] и [к]:

1 -V

Р0п =

=1

Е (1 -V)

(1 + п)(1 - 2п)

(33)

Из первого условия следует V =1/2, что противоречит второму условию. Применим в данном случае приближённый инженерный подход: полагая, что коэффициент Пуассона V близок по величине к 0,5. Точность решения тем выше, чем большее количество девяток позволит ввести программа, реализующая МКЭ. Так, перепишем (33):

IV = 0,499999...

(1 + п)(1 - 2п)

Iе = Р0 п-

(34)

(1 -V)

получим решение уравнения :

юг

81П -

и(г, і) = А-

с

юЬ

СОБ Юі .

(36)

Здесь (о=2пґ — циклическая частота колебаний, Ї — частота колебаний, А — амплитуда возмущения, Ь — длина трубки.

Скорость распространения волны определим по формуле

пЯГ

V м '

где .К = 8,31447 Дж/ (моль'К) — универсальная газовая постоянная, Т — абсолютная температура воздуха в трубке, М = 0,0289644 кг/моль — молярная масса воздуха.

Перепишем амплитуду приращения давления (8) в производных:

л / ч

Ap(z) = -Р0 л —,

Эz

тогда, с учётом (36), можно записать

юг

Лр(г) = -р0 пАю-

юЬ

(37)

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

Это и есть аналитическое решение задачи распространения упругих колебаний в трубке.

Результаты расчётов по МКЭ по [1], по МКЭ в Апзуз и аналитическое решение представлены на (рис. 3, 4).

Результаты расчётов по МКЭ в Апзуз при различной плотности сетки разбиения (20x4 и 40x8) представлены на (рис. 5, 6).

Как видим, результаты расчётов различными методами соответствуют друг другу с достаточной точностью. При этом коэффициент Пуассона V можно

с

V

СОБ

с

с

ОМСКИЙ НАУЧНЫЙ ВЕСТНИК № 1 (107) 2012 МАШИНОСТРОЕНИЕ И МАШИНОВЕДЕНИЕ

МАШИНОСТРОЕНИЕ И МАШИНОВЕДЕНИЕ ОМСКИЙ НАУЧНЫЙ ВЕСТНИК № 1 (107) 2012

ОМСКИЙ НАУЧНЫЙ ВЕСТНИК № 1 (107) 2012 МАШИНОСТРОЕНИЕ И МАШИНОВЕДЕНИЕ

МАШИНОСТРОЕНИЕ И МАШИНОВЕДЕНИЕ ОМСКИЙ НАУЧНЫЙ ВЕСТНИК №1 (107) 2012

-!=200Нг расчёт МКЭ в Апзуэ М=430Н2 расчёт МКЭ в Апэуэ -ЬбООНг расчёт МКЭ в Апэуа -f=600Hz расчёт МКЭ в Апэуэ

"!=200Нг аналитический расчёт по (37) ■1=430Нг аналитический расчёт по (37) ■1=500Нг аналитический расчёт по (37) -1=600Нг аналитический расчёт по (37)

ОМСКИЙ НАУЧНЫЙ ВЕСТНИК № 1 (107) 2012 МАШИНОСТРОЕНИЕ И МАШИНОВЕДЕНИЕ

МАШИНОСТРОЕНИЕ И МАШИНОВЕДЕНИЕ ОМСКИЙ НАУЧНЫЙ ВЕСТНИК № 1 (107) 2012

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

Таким образом, в статье:

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

2. Показана возможность использования изотропного осесимметричного конечного элемента для описания упруго-динамических свойств газа, при этом физические константы элемента определены расчётом по полученным зависимостям (34), что даёт возможность задействовать ресурс программ типа Лмуз и др. для решения задач распространения и взаимодействия с преградами упругих волн.

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

Библиографический список

1. Расчётно-экспериментальные методы проектирования сложных резинокордных конструкций : монография / И. А. Три-бельский [и др.]. — Омск : Изд-во ОмГТУ, 2011. — 240 с.

2. Зенкевич, О. Метод конечных элементов в технике : пер. с англ. / О. Зенкевич. — М. : Мир, 1975. — 544 с.

3. Норри, Д. Введение в метод конечных элементов : пер. с англ. / Д. Норри, Ж. де Фриз. — М. : Мир, 1981. — 304 с.

4. Сегерлинд, Л. Применение метода конечных элементов : пер. с англ. / Л. Сегерлинд. — М. : Мир, 1979. — 392 с.

5. Секулович, М. Метод конечных элементов / М. Секуло-вич : пер. с серб. Ю. Н. Зуева ; под ред. В. Ш. Барбакадзе. — М. : Стройиздат, 1993. — 664 с.

6. Кухлинг, Х. Справочник по физике : пер. с нем. / Х. Кух-линг. — 2-е изд. — М. : Мир, 1985. — 520 с.

ТРИБЕЛЬСКИЙ Иосиф Александрович, доктор технических наук, старший научный сотрудник, заместитель генерального директора по научной работе. БОХАН Владимир Викторович, младший научный сотрудник НПП «Прогресс», аспирант кафедры «Транспорт и хранение нефти и газа, стандартизация и сертификация» Омского государственного технического университета.

Адрес для переписки: [email protected]

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

© И. А. Трибельский, В. В. Бохан

Книжная полка

Трибельский, И. А. Бортовые соединения резинокордных конструкций : монография / И. А. Трибельский; ОмГТУ. - Омск : Изд-во ОмГТУ, 2011. - 130 с. - ISBN 978-5-8149-1109-4.

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

Блюменштейн, В. Ю. Проектирование технологической оснастки : учеб. пособие для вузов по направлению «Конструкторско-технологическое обеспечение машиностроительных производств» / В. Ю. Блюменштейн, А. А. Клепцов. - 2-е изд., испр. и доп. - СПб. [и др.] : Лань, 2011. -219 с. - ISBN 978-5-8114-1099-6.

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

Мухин, В. Ф. Автоматизация сварочных процессов : учеб. пособие для вузов по направлению 150700 «Машиностроение» / В. Ф. Мухин, Е. Н. Ерёмин ; ОмГТУ. - Омск : Изд-во ОмГТУ, 2011. -114 с. - ISBN 978-5-8149-1100-1.

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