CC BY
16
СРАВНИТЕЛЬНАЯ ОЦЕНКА ВАРИАНТОВ ИНТЕРПОЛЯЦИИ ПОЛЕЙ ПЕРЕМЕЩЕНИЙ НА ПРИМЕРЕ ЗАДАЧИ ОСЕСИММЕТРИЧНО НАГРУЖЕННОЙ ОБОЛОЧКИ ВРАЩЕНИЯ
Ю.В. КЛОЧКОВ, д-р техн. наук, профессор А.П. НИКОЛАЕВ, д-р техн. наук, профессор С.С. МАРЧЕНКО, канд. техн. наук, доцент А.А. ШУБОВИЧ, старший преподаватель
Волгоградская государственная сельскохозяйственная академия
г. Волгоград, 400002, пр. Университетский, д. 26;
Тел. (8442) 41-17-84. Факс: (8442) 41-10-85; E-mail: vgsxa@avtlg.ru
Рассмотрены варианты получения интерполяционных выражений на примере конечного элемента осесимметрично нагруженной оболочки вращения с шестью степенями свободы в узле и показана высокая эффективность алгоритма с векторной интерполяцией полей перемещений.
КЛЮЧЕВЫЕ СЛОВА: оболочка вращения, конечный элемент, элемент дискретизации, векторная интерполяция
В настоящее время для анализа напряженно-деформированного состояния большинства оболочечных конструкций строительного и машиностроительного
профиля используют современные вычислительные комплексы такие, как AN-SYS, NASTRAN, ABACUS и другие. Все эти комплексы имеют удобный многофункциональный интерфейс, мощную библиотеку конечных элементов, и в большинстве случаев позволяют эффективно исследовать процессы деформирования конструкций в линейной и нелинейной постановках.
В то же время не следует забывать о том, что адекватность полученных параметров напряженно-деформированного состояния конструкций в определенной степени зависит от тех гипотез, предпосылок, моделей элементов дискретизации, способов аппроксимации и других факторов, которые реализуются в том или ином вычислительном комплексе. Так, например, все современные конечно-элементные расчетные комплексы базируются на элементах, построенных на независимой интерполяционной процедуре, суть которой заключается в том, что каждая компонента вектора перемещения внутренней точки элемента дискретизации интерполируется через узловые значения этой же компоненты и не зависит от узловых значений остальных двух компонент, что не всегда приводит к получению приемлемых результатов.
В настоящей статье рассмотрены варианты получения интерполяционных выражений на примере конечного элемента осесимметрично нагруженной оболочки вращения с шестью степенями свободы в узле, и показана высокая эффективность алгоритма с векторной интерполяцией полей перемещений.
1. Геометрические параметры оболочки
Рассматриваемый элемент дискретизации является фрагментом срединного меридиана с узлами i и j, выделенным двумя плоскостями, перпендикулярными оси вращения (рис. 1). Положение произвольной точки меридиана определяется в декартовой системе координат xoz радиус-вектором
R = xi + r (x) k, (i)
где i , k - орты декартовой системы; x - осевая координата; r (x) - радиус вращения. Базисные векторы произвольной точки срединного меридиана определяются выражениями
ai = R,s = xsi + rxxsk , a = ai x j = -rxxsi + xsk . (2)
С использованием (2) без принципиальных затруднений можно получить соотношения
рт |=иН; а }=нй; {а. }=и{а}, (з)
где {а} = {а1 а} - строка базисных векторов в произвольной точке; {а"} ={аГ а"} - строка базисных векторов узловой точки (т = 1, j);
{а,!} = {а^ а,в} - строка производных базисных векторов произвольной точки;
{а,..} ={а1,!! а,!!} - строка вторых производных базисных векторов.
Глобальные координаты х, . (рис. 1) определяются через узловые значения линейной зависимостью
л = 1] +I]', (4)
2 2 W
где X, Л - дуговые координаты узлов конечного элемента; Л - координата х или !; ] - локальная координата, изменяющаяся в пределах -1 < ] < 1 (рис. 2).
51
Из (4) определяются производные хц , цх, sц и г.
-1
О | г
Рис. 1 Рис. 2
2. Вектор перемещения и его производные
Вектор перемещения произвольной точки срединного меридиана представляется компонентами в базисе этой же точки
: $ {и}.
(5)
v = ua\ + м>а = {и}
1x2 2x1
Дифференцированием (5) с учетом (3) определяются производные
V,, = al (u s + b11u + Ь21^) + a (w , + Ь12и + Ь22= a1tl + at1;
v,ss = al (и,,, + 2u)sbп + иеи + 2w)sb2l + we2l) + (6)
+a (2u,sb12 + ue12 + ^ + 2wsb22 + we22 ) = ^11 + at11.
3. Скалярная аппроксимация перемещений
В качестве узловых неизвестных принимаются векторы - строки в локальной и глобальной системах координат
{чу }г = {ч' ч1 я',г чц ч'гг чгг
{чгу }г = {ч' ч1 ч\* я1 ч'^ ч1^ },
( 7 )
где под символом ч принимается компонента и или w .
Перемещения произвольной точки срединного меридиана определяются интерполяционной зависимостью [1, 2, 3]
и=мТ К}; w={*} {<}, (8)
где {ф}Т - матрица-строка функций формы. Производные компонент вектора перемещения определяются дифференцированием (8)
u,s = ф }Т {иЛ}; w,s = ф }т {<};
u,ss = {ф^}т {<}; ^ = {ф^}г {<}.
Из соотношений (8) и (9) видно, что каждая компонента вектора перемещения и её производные зависят от узловых значений только этой же компоненты.
4. Векторная аппроксимация перемещений
В качестве узловых неизвестных принимаются строки, содержащие векторы перемещений узловых точек и их производные в локальной и глобальной системах координат
\Т /...........4 / „\Т
(9)
~ Л
1 -1
-1 -1
-1
V г V г V,гг V гг
}; (1О) {V; } ={
1 -1 -1 -1 -1 -1
= {V V V, V,s V,ss V,ss
}. (11)
1
На основании дифференциальных зависимостей
v,n = v,ssn ; v^n = V,ss (s,n)
(12)
между векторами (10) и (11) можно сформировать матричное соотношение
}. (13)
Используя (5) и (6), запишем векторы перемещений узловых точек и их производные в виде
-т т - т т - т -т 1ш - т т - т -т 1ш - т т - т
V = ита\ + wma ; V* = /21та1 + t1mа ; V,** = t11mаl + ^а ; (т = i, j ) и сформируем матричное соотношение
;} = [^]Ь } , (14)
6x1
где
A
6x12
—1 a1 a
a1 a
a1 a1
a{ a
a11 a1
a1 —a—j
{ny } = {u w u w tftitij tij tii til t;/t/i}.
1x12
Столбец {ny} с учетом (5), (6) может быть выражен через обычный столбец варьируемых параметров в глобальной системе координат
К} = [P]{<} . (15)
12x1 12x12 12x1
При векторном способе интерполяции перемещений вектор перемещения внутренней точки элемента дискретизации интерполируется через узловые векторы (10), (11)
v=мТ {vy }=мТ {v;}=wr як}, (16)
где элементы строки {/}Т определены с использованием (13) выражениями
/1 = ^; / = ^; гз = V3s'n; г4 = пs,n; У5 = % (s,n) ; г6 = % (s,n) .
Производные вектора перемещения (16) определяются выражениями v,s = Ь}Т [L][A]{ny} = {/,}Т [A]{ny}; V,ss = {<P,ss }Т [ L][A]K } = {/,„ }Т [ A]{ny }. После замены базисных векторов узловых точек в матрице [ а] через базисные векторы внутренней точки конечного элемента по (3) соотношения (16) и (17) можно представить в виде
(17 )
v = {«}Т [ H ] {ny}; v,s = {«}Т [ Hs ] {ny}; v,ss = {я}Т [ H ss ] {ny } :
1x2 2x12 12x1
1x2 2x12 12x1
1x2 2x12 12x1
(18)
6x12
где [Н]=
2x12
71 [<? ] У2 ] Уз [<? У4 У5 [>' Уб ]7
[ Н, ] =
2x12
[ Н.к ] =
У1,
М
2x2
2x2 Т
У2,
■[л11
2x2
2x2 Т
Уз,
■[л" ]Т
2x2
У4,
]Т
2x2
У5,
:[л']Т
2x2
Уб,
^ ]Т
2x2
;(19)
У1,
,[✓ ]Т
2x2
У2,
■[л1 ]Т
2x2
Уз,
:[л']Т
2x2
У4,
^ ]Т
2x2
У5,;
:[л']Т
2x2
Уб,
■[л1 ]Т
2x2
Приравнивая правые части (5), (б) и (18), получим интерполяционные выражения
и=ИйЛ; 1? }=Н 1К}; }?' }=[н,- 1К}, (20)
2x12 12 x1 L 1 J 2x12 12x1 L 11J 2x12 12x1
которые используются при получении матрицы жесткости конечного элемента. Анализируя интерполяционные зависимости (20), можно отметить, что при
векторном способе интерполяции каждая из компонент вектора перемещения внутренней точки конечного элемента зависит от узловых значений как вектора {иГ} ,
так и вектора {^Г } .
Пример расчета. В качестве примера рассмотрим напряженно-^ деформированное состояние осе-симметричной оболочки вращения, радиус которой задавался функциональной зависимостью вида
(рис. 3)
Рис. 3
г = А + В ^ (х/С), А = 1,3м; В = 0,4м .
Таблица 1
2x2
2x2
2x2
2x2
2x12
Координата х, м ^^ Сетка а , МПа^\ 1x17 1 x 33 1 x 49 1 x 65
0,00 12,228 12,214 12,234 12,233 12,235 12,234 12,235 12,235
авв -3,706 -3,710 -3,793 -3,793 -3,808 -3,808 -3,814 -3,814
0,48ж 0,009 -0,002 0,001 0,000 0,000 0,000 0,000 0,000
авв 17,534 17,531 17,531 17,531 17,530 17,530 17,530 17,530
Оболочка была нагружена внутренним давлением интенсивности д = 0,2МПа ; толщина оболочки ; = 0,01 м; модуль упругости Е = 2,06-105МПа; коэффициент Пуассона у = 0,3. Левый край оболочки был шарнирно закреплен, а правый край был свободен.
Первоначально величина параметра С была принята равной С = 0,48 м , а осевая координата x изменялась в пределах 0 < х < 0,48ж м .
Таблица 2
Координата X , м \Сетка \ уз- \лов а , \ МПа\ 20 X 20 30 x 30 40 x 40 60 x 60
shell 63 shell 93 shell 63 shell 93 shell 63 shell 93 shell 63 shell 93
0,00 12,19 12,03 12,47 11,84 12,17 12,11 12,37 11,95 12,15 12,15 12,29 12,04 12,14 12,18 12,21 12,12
-3,74 -3,80 -3,72 -3,91 -3,77 -3,81 -3,76 -3,84 -3,78 -3,79 -3,77 -3,84 -3,78 -3,79 -3,78 -3,83
0,48п &XX 0,064 1,32 0,034 -0,03 0,034 0,092 0,036 -0,03 0,023 0,068 0,027 -0,03 0,00 0,059 0,013 -0,01
17,49 17,63 17,54 17,17 17,50 17,54 17,49 17,14 17,49 17,51 17,45 17,12 17,48 17,47 17,41 17,10
При данном значении параметра C кривизна оболочки в концевых сечениях равна
K = |г)ХСх3х| = 5,64 м-,
что позволяет говорить о том, что рассчитываемая оболочка достаточно пологая. Расчеты были выполнены в двух вариантах. В первом варианте был реализован описанный выше алгоритм, основанный на векторном способе интерполяции перемещений. В качестве контрольного варианта применялся программный комплекс ANSYS, причем в последнем для решения поставленной задачи использовались два типа элемента: «shell 63 линейный» и «shell 93 параболический». Результаты расчета представлены в форме таблиц, в которых приведены значения меридиональных и кольцевых напряжений на внутренней и наружной поверхностях оболочки в концевых сечениях в зависимости от сетки дискретизации оболочки. Табл. 1 соответствует первому варианту расчета, а табл. 2 -контрольному.
Сопоставительный анализ числовых данных, представленных в табл. 1, 2 показывает, что при расчете пологой оболочки параметры напряженно- деформируемого состояния в обоих вариантах практически одинаковы. Кроме того, следует отметить, что численные значения напряжений в контрольном варианте ANSYS, полученные при использовании элементов «shell 63 линейный» и «shell 93 параболический» также совпадают между собой.
Для верификации полученных значений напряжений можно использовать условие равновесия, согласно которому меридиональное напряжение в точке шарнирного опирания может быть вычислено по формуле
2 2 пп - пг2 .
= —^—- • qfr,
2nrj
где г и г2 - радиусы вращения в точке шарнирного опирания (х = 0) и в концевом сечении (х = 0,48ж м ) оболочки соответственно. При подстановке исходных данных получим
л .1 72 _л.о о2 ам =------0,2/0,01 = 12,235МПа .
2 л ■ 1,7
Как видно, полученное из условия равновесия значение меридионального напряжения практически совпадает с численными значениями, полученными в процессе конечно-элементного решения в обоих вариантах расчета. Кроме того, следует отметить устойчивую сходимость вычислительного процесса при увеличении числа элементов дискретизации.
Если параметр С уменьшить в 6 раз и принять равным С = 0,08 м , а значение осевой координаты х изменять в пределах 0 < х < 0,08л м, то кривизна оболочки значительно возрастет, и её нельзя уже будет рассматривать как пологую. Однако, расчетная схема оболочки (см. ис. 3) останется прежней и для неё также будет справедливо рассмотренное выше условие равновесия.
Таблица 3
Координата х, м Сетка ^-улов а , МПа"\ 1 х 81 1 х 101 1 х 161 1 х 201
0,00 13,59 9,37 12,94 10,73 12,41 11,86 12,32 12,04
авв -700,3 -701,6 -700,6 -701,2 -700,8 -701,0 -700,8 -700,9
0,08л 0,392 -0,186 0,206 -0,097 0,052 -0,024 0,02 -0,012
авв -2,452 -2,625 -2,494 -2,584 -2,549 -2,572 -2,623 -2,635
Таблица 4
Координата х, м \ Сетка \узлов а, \ МПа \ 90 х 60 150 х 50 200 х 50
shell 63 shell 93 shell 63 shell 93 shell 63 shell 93
0,00 ахх 26,13 20,97 8,99 3,78 22,65 19,58 9,47 7,83 21,40 20,42 8,58 8,92
-757,4 -756,7 -642,1 -933,5 -752,4 -753,4 -654,6 -915,8 -752,7 -752,6 -654,6 -916,1
0,08л ахх 2,31 1,37 1,38 1,37 1,77 0,69 0,07 0,79 1,89 0,59 -0,22 0,52
5,19 4,60 -5,17 -9,15 4,66 4,63 -7,14 -10,98 4,49 4,55 -8,4 -12,00
Результаты конечно-элементных расчетов при уменьшенном значении параметра С приведены в таблицах 3 и 4, структура которых аналогична табл. 1 и 2. Сопоставительный анализ численных значений напряжений, приведенных в табл. 3 и 4 показывает, что контролируемые параметры напряженно- деформируемого состояния при измененном значении параметра С существенно различаются между собой.
Так, например, численные значения меридионального напряжения в точке шарнирного опирания в первом варианте монотонно приближаются к значению, вычисленному из условия равновесия с увеличением числа элементов дискретизации. Кроме того, следует отметить достаточно быструю сходимость значений кольцевых напряжений в рассматриваемых сечениях.
В контрольном варианте ANSYS следует отметить значительное различие в численных значениях напряжений, полученных при использовании элементов «shell 63» и «shell 93», чего не наблюдалось при расчете пологой оболочки. Так, например, меридиональные напряжения в точке шарнирного опирания различаются между собой более чем в два раза и не совпадают со значениями aM = 12,235МПа, полученные из условия равновесия. Кроме того, окружные напряжения на правом краю оболочки имеют разные знаки.
Таким образом, можно сделать вывод, что при расчете непологих оболочек программный комплекс ANSYS не позволяет получать адекватную картину напряженно-деформируемого состояния.
При расчете непологих оболочек следует использовать элементы дискретизации, матрицы жесткости которых сформированы на основе векторного способа интерполяции перемещений.
Л и т е р а т у р а
1. Зенкевич О.М. Метод конечных элементов в технике. / О.М. Зенкевич. - М.: Мир.-1975.- 542 с.
2. Постное В.А., Хархурим И.Я. Метод конечных элементов в расчетах судовых конструкций. Л.: Судостроение, 1974. 344 с.
3. Голованов А.И. Современные конечно-элементные модели и методы исследования тонкостенных конструкций / А.И. Голованов, А.В. Песошин, О.Н. Тюленева.- Казань: Изд-во КГУ.- 2005.- 442 с.
4. Николаев А.П. Расчет оболочек на основе МКЭ в двумерной постановке / Николаев А.П., Клочков Ю.В., Киселев А.П.,. Гуреева Н.А.- Волгоград: ИПК ФГОУ ВПО ВГСХА «Нива».- 2009.- 196 с.
COMPARATIVE APPRASAL OF VARIANTS OF INTERPOLATION DISPLACEMENT FIELDS ON THE EXAMPLES OF THE PROBLEM AXIALLY SYMMETRIC THE LOADED SHELL OF REVOLUTION
J.V. Klochkov, A.P. Nikolaev, S.S. Marchenko, A.A. Shubovich
Now for the analysis mode of deformation of the majority shell of constructions of a building and machine-building structure use modern computer complexes such, as ANSYS, NASTRAN, ABACUS and others.
All these complexes have the convenient multipurpose interface, powerful library of finites elements, and in most cases effectively allow to investigate processes of deformation of constructions in linear and nonlinear statements.
At the same time it is not necessary to forget that adequacy of the received parameters is mode of deformation of constructions in the certain degree depends on those hypotheses, preconditions, models of elements of discretization, ways of approximation and other factors which are realized in this or that computer complex.
So, for example, all modern finites elements settlement complexes are based on the elements constructed on independent interpolation to procedure which essence consists that everyone of a component of a vector of displacement of an interior point of an element of discretization is interpolated through pivotal values of this components and other two component does not depend on modes values, that not always reduce to reception of acceptable results.
In present article view variants of reception interpolation expressions on an example of a finite element axially symmetric the loaded shell of rotation with six degrees of freedom in node and are considered, and high efficiency of algorithm with vector interpolation of displacement fields.
KEY WORDS: shell of revolution, finite element, vector interpolation.
