CC BY
147
Б01: 10.15593/2224-9982/2015.40.06 УДК 539.3: 629.7.023
Р.В. Бульбович, В.В. Павлоградский, П.П. Еременко
Пермский национальный исследовательский политехнический университет, Пермь, Россия
РАСЧЕТ НАПРЯЖЕННО-ДЕФОРМИРОВАННОГО СОСТОЯНИЯ ТЕЛА ИЗ УПРУГОГО ОРТОТРОПНОГО МАТЕРИАЛА С ПОМОЩЬЮ МЕТОДА КОНЕЧНЫХ ЭЛЕМЕНТОВ
Рассматривается расчет напряженно-деформированного состояния (НДС) линейного упругого ортотропного материала с привлечением метода конечных элементов, реализованный при помощи языка программирования С++. Используется подход «препроцессор ДЫЗУЗ - собственный решатель - постпроцессор ДЫЗУЗ». Показана возможность использования ДЫЗУЗ для автоматизации формирования конечно-элементной модели при помощи языка параметрического программирования ДРйЦ а также визуализации и анализа результатов расчета собственной программы при помощи постпроцессора ДЫЗУЗ. Приводятся физические соотношения для линейного упругого ортотропного материала в системе координат, совпадающей с осями симметрии упругих свойств материала. Рассмотрены особенности задания коэффициентов Пуассона для ортотропного материла в зависимости от их положения в записи обобщенного закона Гука. Выводится общий вид матрицы упругости ортотропного материала в осях симметрии упругих свойств для общего случая трехмерного тела. Особое внимание уделяется алгоритму пересчета матрицы упругости материала в систему координат, не совпадающую с осями симметрии упругих свойств. В качестве примера приводятся результаты расчета НДС конической оболочки, моделирующей в первом приближении раструб соплового блока ракетного двигателя на твердом топливе, оси симметрии упругих свойств которой расположены под углом к оси вращения оболочки. Результаты расчета сравниваются с результатами, полученными в пакете ДЫЗУЗ. Для вычисления значений деформаций и напряжений в узлах конечно-элементной модели используется теория сопряженной аппроксимации. Проведено сравнение результатов расчета в собственной программе и в пакете ДЫЗУЗ по значениям узловых перемещений, деформаций и напряжений. При расчете узловых перемещений получена хорошая сходимость с расчетом в ДЫЗУЗ, при расчете узловых деформаций и напряжений сходимость хуже, что свидетельствует о необходимости дальнейшего тестирования метода сопряженной аппроксимации и поиска и анализа других методов расчета узловых деформаций и напряжений.
Ключевые слова: ракетный двигатель твердого топлива, метод конечных элементов, анизотропный материал, численное моделирование, напряженно-деформированное состояние, ДЫЗУЗ.
R.V. Bulbovich, V.V. Pavlogradskiy, P.P. Eremenko
Perm National Research Polytechnic University, Perm, Russian Federation
THE CALCULATION OF THE STRESS-STRAIN STATE OF ELASTIC ORTHOTROPIC BODY BY FINITE ELEMENT METHOD
The article discusses calculation of stress-strain state of a linear elastic orthotropic material by the finite element method using the C++ programming code. The approach "ANSYS preprocessor -own solver - ANSYS postprocessor" is used. It is shown the ability to use ANSYS for automatic formation of finite element models using parametric programming language APDL, as well as using the ANSYS postprocessor for visualization and analysis of results calculated by own programs. Physical relations are presented for linear elastic orthotropic material in the coordinate system coinciding with the symmetry axes of the elastic properties of the material. The peculiarities of assignment of the Poisson ratios for an orthotropic material depending on their position in generalized law of Hooke are considered. The matrix of elastic orthotropic material in symmetry axes of the elastic properties for the general case of a three-dimensional body is derived. Special attention is paid to the recalculation algorithm of the elasticity matrix for the coordinate system which does not coincide with the symmetry axes of the elastic properties. As an example the calculated stress-strain state of a conical shell simulating in a first approximation a solid-propellant rocket motor nozzle is presented. The symmetry axis of the elastic properties are set at an angle to the axis of rotation of the shell. The calculated results are compared with those obtained in the package ANSYS. To calculate the values of strains and stresses at the nodes of the finite element model the theory of conjugate approximation is used. Nodal values of displacements, strains and stresses obtained by own program are compared with results obtained by ANSYS. The nodal displacements are found to be in a good agreement, but accuracy of strains and stresses is worse. It evidences that it is necessary further testing the conjugate approximation method and search other methods of calculation of nodal strains and stresses.
Keywords: solid-propellant rocket engine, finite element method, anisotropic material, numerical modeling, stress-strain state, ANSYS.
Создание сложных конструкций в ракетной технике стало возможным прежде всего благодаря тщательному анализу их элементов, базирующемуся на современных методах расчета. В настоящее время для расчета напряженно-деформированного состояния (НДС) элементов конструкции очень эффективным оказывается метод конечных элементов (МКЭ).
Конструкция современных ракетных двигателей твердого топлива (РДТТ), в частности сопловых блоков, включает как металлические элементы, так и значительную часть элементов из композиционных материалов. При оценке НДС конструкции при работе РДТТ необходимо учитывать пластические свойства материалов (как металлов, так и композитов), в противном случае расчет дает завышенные значения напряжений, при которых элементы конструкции разрушаются. Этого не происходит именно за счет пластических свойств материалов.
Актуальной проблемой является моделирование пластических свойств композитов, которые отличаются от металлов тем, что проявляют анизотропию свойств, т.е. каждому направлению соответствует своя диаграмма деформирования, которая, кроме того, зависит от компонент напряженного состояния.
В статье рассматривается один из вопросов, возникающих при решении данной задачи, а именно расчет НДС для ортотропного линейно-упругого материала. Результаты этого расчета будут использованы в дальнейшем для моделирования пластического поведения композитов.
Расчет проводится для конической трубы переменной толщины под действием внутреннего давления, меняющегося по заданному закону (в данном случае линейному) по длине образующей. Задача решается в осесимметричной постановке при статическом нагружении конструкции. Размеры расчетной области приведены на рис. 1. Здесь гв = 50 мм, И1 = 30 мм, И2 = 10 мм, а = 15°.
В расчете используется четырехугольный восьмиузловой субпараметрический конечный элемент (элемент имеет прямолинейные границы). При расчете НДС учтена ориентация осей симметрии ортотропного материала. Принимается, что оси симметрии материала направлены согласно рис. 2, где т1Х111 - система координат, совпадающая с осями симметрии ортотропного материала (в данном случае ось г1 параллельна образующей внутренней конической поверхности). Конечно-элементная сетка, положение осей симметрии материала и граничные условия показаны на рис. 2.
\
1 а
1
Рис. 1. Размеры расчетной области
г
Рис. 2. Конечно-элементная сетка, положение осей симметрии материала и граничные условия
Данные для формирования конечно-элементной модели (конечно-элементная сетка, номера узлов, в которых запрещены перемещения, номера элементов, к которым прикладывается давление) формируются автоматически в пакете АКБУБ при помощи макроса на языке АРБЬ [1-3], создающего текстовые файлы с исходными данными.
Упругие свойства материала, заданные в расчете, приведены ниже.
Упругие свойства материала:
где Ег, Еи Ег - модули упругости соответственно в направлениях г1, г1, Ь; Оп, Ой, Огг - модули сдвига соответственно в плоскостях г^, Х111, г1г1; цгг, цгг - коэффициенты Пуассона ([М/ - абсолютная
величина отношения относительного удлинения в направлении г к относительному поперечному сжатию в направлении у). Избыточные характеристики выделены жирным шрифтом.
Ортотропный материал имеет девять независимых упругих постоянных. Между упругими постоянными имеется связь:
Ег = 5103 МПа, Е = 200 1 03 МПа, Ег = 50 103 МПа, Оп = 12,5103 МПа, О ь = 25 103 МПа, ОГ1 = 50 103 МПа,
Цгг = 0,1,
Мй = 0,15, М г = 0,0075, Ръ = 0,01, Рп = 0,6,
= 0,0001875.
В связи с тем, что при расчете матриц жесткости элементов и, как следствие, матрицы жесткости конструкции используются координаты узлов в системе координат г^, а свойства материала заданы в системе координат Т\1\2\, возникает необходимость перевода матрицы упругости элемента [Б] в систему координат Н2. Для этого используются следующие формулы [4]:
- £ б ъП
= ЩЩ £ £ -ЧтЯп. ,
ш=! п=! ЮтЮ„
где Б. - матрица упругости размером 6*6 для общего случая трехмерного тела в системе координат г^; Опп - матрица упругости размером 6*6 для общего случая трехмерного тела в системе координат г!1!1!; юк = ! при к = 4, 5, 6, юк = 2 при к = !, 2, 3; я - матрица размером
6*6, имеющая вид
(
Ч
а2 а2 а2 2а2а3 2аза1 2а!а2
Р2 Р2 в? 2р2рз 2рзр! 2р!р2
У2 у2 У? 2У2 Уз 2Уз У! 2У!У2
Рл в2 У 2 Рз Уз р2Уз + рзУ2 Р!Уз + рзУ2 Р!У2 + Р2У!
У!а! У2а2 Узаз У2аз + Уза2 У1аз + Уза2 У!а2 + У2а1
«А а2в2 азрз а2рз + «зр2 а!рз + азр2 а!р2 + а2р1 )
где аг, Рг, Yi - направляющие косинусы осей системы координат г!^!г! в системе координат
Приведенные выше формулы применимы для общего случая трехмерного тела. Для того чтобы воспользоваться ими в случае осе-симметричной задачи, необходимо сначала получить матрицу [ Б ] для трехмерного случая, используя заданные упругие характеристики, а затем удалить из нее строки и столбцы, соответствующие касательным напряжениям и тгг, получив при этом матрицу размером 4*4.
Для получения выражения матрицы упругости через характеристики материала запишем обобщенный закон Гука для ортотропного материала в осях симметрии [5, 6]:
1 а Мг а М2г а
- °г----С 7 ,
Е " Е, ' Е7 7
в = М г, а + 1 а М 7, а в, ---аг +--а,--а7,
' Ег г Е ' Е„ 7
1
в7 =-Ь7аг - ^а, + —аг,
7 Ег Е, Ег
1,7
О7
=
У г7 ✓—Т Тг7 ,
Ог7
=_
У г, = ^ Тг,.
Ог,
Выражая напряжения через деформации и используя обозначения (Вг в, в2 уУг7 Уг,)Т и {а дующее выражение для [ й ]:
{в} = (Вг в, в2 у Уг7 у г,)Т и {а} = (аг а, а2 т ^ т^ тп)т [7-9], получим сле
й --I
( Ег (("М* -1) -Ег (Мг + М7г М,7 Ег (М7г + М,г "М* ) 0
( + М г7 М7, ) Е (Мг7 ' М 7г -1) ( +Мг, * М 7г ) 0
-Е7 ( + Мг, М,7 )-Е7 (( + Мг -Мг7 ) Е7 ( -Мг -1) 0
0 0 0 АО
0 0 0
0 0 0 0
0 0 0 0
7
0 АО,,,
0 0 0 0 0
0 А-О
гг /
где А = Мг,• М,г + М7 Мгг + М7 М7, + Мг,• Мгг' М,7 + М„' Мгг' М7, - 1.
Тогда окончательный вид матрицы упругости для осесиммет-ричного случая в осях симметрии материала запишется следующим образом:
' Ег (,7 М7, -1) -Ег (,г + М7г М,7 ) -Ег (М7г + М,г ' М7, ) 0 1 -Е, (г, + Мг7 М7, ) Е (Мг7 " М7г -1) -Е, ( + Мг, ' М7г ) 0 А -Е7 (Мг7 + Мг, " М,7 ) -Е7 (М,7 + М,г " Мг7 ) Е7 (Мг, ' М,г -1) 0
0 0 0 А-ОГ7
При задании нагрузки, меняющейся вдоль внутренней конической поверхности, принимается, что на каждый конечный элемент действует нагрузка, заданная линейным законом вдоль его границы. Значения узловых сил вычислялись по формулам
ловых сил, действующих в узлах нагруженной поверхности конечного элемента (номерам 1 и 3 соответствуют крайние узлы, номеру 2 - узел посередине); г1, 71, г3, 73 - координаты крайних узлов; р1, р3 - величина давления в крайних узлах.
Решение задачи реализовано в виде программы на языке программирования С++.
Результаты расчета получены в виде текстовых файлов с массивами, содержащими перемещения узлов, напряжения и деформации, приведенные к узлам. Для расчета значений деформаций и напряжений в узлах использовалась теория сопряженной аппроксимации [10].
30
где /}, /], /у, /7, /г3, /7 - радиальные и осевые составляющие уз-
Для визуализации результатов расчета программы использовался постпроцессор пакета АКБУБ. На рис. з показаны поля радиальных и осевых перемещений на деформированной форме тела; на рис. 4 - напряжения в теле. Значения перемещений и напряжений показаны в системе координат
.008593 .075703 .142814 .209924
.042148 .109259 .176369 .260257
а
.368919 -.270541 -.172162 -.073784
-.31973 -.221351 -.122973 .555Е-16
б
Рис. з. Перемещения: а - радиальные; б - осевые
-173.916 126.779 427.475 728.17
-23.5683 277.127 577.823 953.692
-400.876 -292.053 -183.231 -74.408
-346.465 -237.642 -128.819 7.209
-30.993 -15.8231 -.653267 14.5166
-23.4081 -8.2382 6.93167 25.894
Рис. 4. Напряжения: а - радиальные; б - окружные; в - осевые; г - касательные
а
б
в
г
Для сравнения результатов расчета в разработанной программе с результатами расчета в пакете АКБУБ приведем значения узловых перемещений, деформаций и напряжений в узлах с номерами от 1 до 10 (всего в модели 1149 узлов). Величины узловых перемещений приведены в табл. 1, деформаций - в табл. 2, напряжений - в табл. 3.
Таблица 1
Перемещения в узлах конечных элементов, мм
Номер Разработанная программа АШУЯ
узла иг и иг иг
1 0,248471 0,000000 0,248527 0,000000
2 0,008593 -0,358851 0,008587 -0,358863
3 0,254264 -0,011335 0,254329 -0,011364
4 0,257571 -0,021664 0,257597 -0,021700
5 0,259258 -0,030543 0,259227 -0,030554
6 0,260068 -0,038647 0,260045 -0,038611
7 0,260345 -0,046043 0,260344 -0,046040
8 0,260257 -0,052925 0,260249 -0,052921
9 0,259934 -0,059383 0,259927 -0,059383
10 0,259447 -0,065512 0,259447 -0,065510
Таблица 2
Деформации в узлах конечных элементов
Номер узла Разработанная программа АШУЯ
£г е* ег е* ег
1 -0,0270 0,0050 -0,0071 0,0114 -0,0272 0,0049 -0,0070 0,0116
2 0,0000 0,0001 0,0000 0,0000 0,0000 0,0001 0,0000 0,0000
3 -0,0275 0,0050 -0,0069 0,0118 -0,0274 0,0050 -0,0068 0,0118
4 -0,0274 0,0051 -0,0066 0,0119 -0,0276 0,0050 0,0000 0,0121
5 -0,0271 0,0050 -0,0062 0,0121 -0,0272 0,0050 -0,0062 0,0121
6 -0,0269 0,0050 -0,0058 0,0120 -0,0268 0,0050 -0,0058 0,0121
7 -0,0266 0,0050 -0,0056 0,0121 -0,0266 0,0049 -0,0056 0,0121
8 -0,0263 0,0049 -0,0054 0,0120 -0,0263 0,0049 -0,0054 0,0120
9 -0,0261 0,0049 -0,0052 0,0120 -0,0260 0,0048 -0,0052 0,0120
10 -0,0259 0,0048 -0,0050 0,0120 -0,0258 0,0048 -0,0050 0,0120
Таблица 3
Напряжения в узлах конечных элементов, МПа
Номер узла Разработанная программа АШУЯ
ог о* О Ог о* ог
1 -146,14 896,05 -158,49 -6,06 -147,11 888,98 -151,44 -4,51
2 -0,27 16,63 0,28 -0,35 0,13 16,59 0,14 0,03
3 -147,54 919,55 -144,02 -1,76 -146,84 908,47 -137,13 -0,26
Окончание табл. 3
Номер узла Разработанная программа А№У8
Сг С, сг 'Тг2 Сг С, Сг
4 -146,67 935,61 -124,68 1,49 -146,57 927,97 -122,83 3,99
5 -141,96 947,39 -104,60 10,02 -143,16 935,10 -105,27 8,70
6 -141,72 953,69 -83,62 10,31 -139,76 942,23 -87,72 13,42
7 -137,65 953,69 -73,51 16,12 -137,36 939,73 -76,67 16,18
8 -136,23 951,32 -62,67 17,25 -134,96 937,23 -65,61 18,94
9 -133,36 946,50 -56,47 21,00 -133,10 931,41 -58,86 20,99
10 -131,95 941,12 -49,15 22,45 -131,24 925,60 -52,10 23,04
Анализ результатов расчета разработанной программы и пакета АКБУБ показал отличную сходимость результатов по узловым перемещениям (различие не более 0,3 %), что говорит о правильности расчета упругих характеристик в системе координат, отличной от осей симметрии свойств материала, и в целом о высокой точности решения задачи.
Расчет узловых деформаций показал удовлетворительную сходимость. Для большинства узлов различие составило не более 2 %, однако в нескольких узлах модели различие достигло ~ 100-200 %, что в целом не повлияло на картину распределения деформаций ввиду малости абсолютных величин деформаций в этих узлах по сравнению с окружающими.
Расчет узловых напряжений также показал удовлетворительную сходимость. Для большинства узлов различие составило не более 13 %, однако в нескольких узлах модели оно достигло ~ 100 %, что также в целом не повлияло на картину распределения напряжений ввиду малости абсолютных величин напряжений в этих узлах по сравнению с окружающими.
Следует отметить, что узловые перемещения определены с высокой точностью, следовательно, недостаточно высокая точность определения узловых деформаций и напряжений обусловлена выбором метода их расчета, т.е. методом сопряженной аппроксимации. В инженерной практике при расчете на прочность в большинстве случаев искомыми величинами являются именно напряжения и деформации в конструкции. В связи с этим актуальным остается вопрос выбора метода расчета напряжений и деформаций через узловые перемещения, дающего более точные значения, конечно, при условии, что перемещения определены верно.
При задании коэффициентов Пуассона ортотропного материала следует обращать внимание на различие коэффициентов и j а также на их величины, которые для ортотропного материала могут быть больше 0,5. Данные коэффициенты должны быть заданы таким образом, чтобы матрица упругости была положительно определена.
Рассмотренный выше метод расчета упругих характеристик в системе координат, отличной от осей симметрии свойств материала, может быть использован и в общем случае анизотропного тела (т.е. с учетом коэффициентов взаимного влияния и коэффициентов Чен-цова), что может быть полезным, например, при расчете НДС элементов РДТТ сложной формы, изготовленных из композиционных материалов.
Выражения для вычисления узловых сил в элементе, нагруженном линейно изменяющимся вдоль его стороны давлением, в дальнейшем будут использованы для задания газодинамических нагрузок, меняющихся вдоль соплового блока. При этом формулы являются универсальными, т.е. позволяют рассчитывать узловые силы для элемента, ориентированного произвольным образом.
Разработанная программа и полученные результаты будут в дальнейшем использованы для расчета НДС тел сложной формы, проявляющих криволинейную анизотропию (раструбы сопловых блоков, корпуса и пр.). Кроме того, расчет НДС упругого ортотропного материала будет использован в итерационной процедуре расчета НДС ортотропного упругопластического материала (при помощи метода переменных параметров упругости).
Библиографический список
1. ANSYS Parametric Design Language Guide. ANSYS Release 14.5. ANSYS Inc., 2012.
2. Каплун А.Б., Морозов Е.М. Олферьева М.А. ANSYS в руках инженера: практическое руководство. - М.: Едиториал УРСС, 2003. - 272 с.
3. Басов К.А. ANSYS: справочник пользователя. - М.: ДМК Пресс, 2005. - 640 с.
4. Лехницкий С.Г. Теория упругости анизотропного тела. - М.; Л.: Гос. изд-во технико-теор. лит-ры, 1950. - 299 с.
5. Ташкинов А.А. Упругость анизотропных материалов: конспект лекций. - Пермь: Изд-во Перм. нац. исслед. политехн. ун-та, 2010. - 49 с.
6. Темопрочность деталей машин / под ред. И.А. Биргера, Б.Ф. Шорра. - М.: Машиностроение, 1975. - 455 с.
7. Зенкевич О. Метод конечных элементов в технике. - М.: Мир, 1975. - 560 с.
8. Усюкин В.И. Строительная механика конструкций космической техники. - М.: Машиностроение, 1988. - 392 с.
9. Бате К., Вилсон Е. Численные методы анализа и метод конечных элементов. - М.: Стройиздат, 1982. - 448 с.
10. Сегерлинд Л. Применение метода конечных элементов. - М.: Мир, 1979. - 392 с.
References
1. ANSYS Parametric Design Language Guide. ANSYS Release 14.5. ANSYS Inc., 2012.
2. Kaplun A.B., Morozov E.M., Olfereva M.A. ANSYS v rukakh in-zhenera: prakticheskoe rukovodstvo [ANSYS in engineer's hand]. Moscow: Editorial URSS, 2003. 272 p.
3. Basov K.A. ANSYS: spravochnik polzovatelya [ANSYS: user directory]. Moscow: DMK Press, 2005. 640 p.
4. Lekhnitskiy S.G. Teoriya uprugosti anizotropnogo tela [Theory of anisotropic elasticity]. Moscow, Leningrad: Gosudarstvennoe izdatelstvo tekhniko-teoreticheskoy literatury, 1950. 299 p.
5. Tashkinov A.A. Uprugost anizotropnykh materialov [Elasticity of anisotropic materials]. Permskiy natsionalnyy issledovatelskiy politekhni-cheskiy universitet, 2010. 49 p.
6. Birger I.A., Shorr B.F. Termoprochnost detaley machin [Thermal strength of machine elements]. Moscow: Mashinostroenie, 1975. 455 p.
7. Zenkevich O. Metod konechnykh elementov v tekhnike [Finite element method in technology]. Moscow: Mir, 1975. 560 p.
8. Usyukin V.I. Stroitelnaya mekhanika konstruktsiy kosmicheskoy tekhniki [Structural mechanics of space technology constructions]. Moscow: Mashinostroenie, 1988. 392 p.
9. Bathe K., Wilson E. Chislennye metody analiza i metod konechnykh elementov [Numerical methods in finite element analysis]. Moscow: Stroyizdat, 1982. 448 p.
10. Segerlind L. Primenenie metoda konechnykh elementov [Applied finite element analysis]. Moscow: Mir, 1979. 392 p.
Об авторах
Бульбович Роман Васильевич (Пермь, Россия) - доктор технических наук, профессор кафедры «Ракетно-космическая техника и энергетические системы» ФГБОУ ВПО ПНИПУ (614990, г. Пермь, Комсомольский пр., д. 29, e-mail: dekan_akf@pstu.ru).
Павлоградский Виктор Васильевич (Пермь, Россия) - кандидат технических наук, доцент кафедры «Ракетно-космическая техника и энергетические системы» ФГБОУ ВПО ПНИПУ (614990, г. Пермь, Комсомольский пр., д. 29, e-mail: rkt@pstu.ru).
Еременко Петр Петрович (Пермь, Россия) - аспирант, ассистент кафедры «Ракетно-космическая техника и энергетические системы» ФГБОУ ВПО ПНИПУ (614990, г. Пермь, Комсомольский пр., д. 29, e-mail: ka-ramba@mail.ru).
About the authors
Roman V. Bulbovich (Perm, Russian Federation) - Doctor of Technical Sciences, Professor, Department of Rocket and Space Engineering and Power Generating Systems, Perm National Research Polytechnic University (29, Komsomolsky av., Perm, 614990, Russian Federation, e-mail: dekan_akf@pstu.ru).
Victor V. Pavlogradskiy (Perm, Russian Federation) - Ph. D. in Technical Sciences, Associate Professor, Department of Rocket and Space Engineering and Power Generating Systems, Perm National Research Polytechnic University (29, Komsomolsky av., Perm, 614990, Russian Federation, e-mail: rkt@pstu.ru).
Petr P. Eremenko (Perm, Russian Federation) - Doctoral Student, Assistant Lecturer, Department of Rocket and Space Engineering and Power Generating Systems, Perm National Research Polytechnic University (29, Komsomolsky av., Perm, 614990, Russian Federation, e-mail: ka-ramb a@mail.ru).
Получено 16.01.2015
