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

Полный текст

Введение Теория механики сплошной среды и, в частности, оболочек разработана достаточно полно [1-3]. Уравнения для определения напряженно-деформированного состояния конструктивных элементов агропромышленного комплекса, химического и авиационного машиностроения и других отраслей получились по сложности такими, что их использование в практике инженерных расчетов оказалось весьма ограниченным. Из-за сложности получения аналитических решений значительное количество исследований посвящено разработке численных и приближенных методов расчета деформируемых тел [4-8]. Среди методов численного определения напряженно-деформированного состояния (НДС) инженерных структур различного назначения метод конечных элементов (МКЭ) является одним из наиболее эффективных, что отражено в работах как отечественных [9-14], так и зарубежных исследователей[21] [15-22]. Определение НДС тонкостенных конструкций в указанных работах выполнялось при использовании гипотезы прямой нормали на основе МКЭ в формулировке метода перемещений. Причем в криволинейных системах координат традиционные функции формы использовались непосредственно для аппроксимации компонент вектора перемещений через узловые величины, в результате каждая компонента вектора перемещений внутренней точки конечного элемента выражалась только через узловые значения этой же компоненты[22] [15; 17-22] и не зависела от других компонент. Эти положения использовались в матрицах жесткостей конечных элементов и при исследовании процессов деформирования в геометрически нелинейной постановке [19; 20; 22; 23]. На основе МКЭ созданы и широко применяются универсальные коммерческие программные продукты типа ANSYS, NFSTRAN, LS-DYNA, ADINA, ASTRA-MOBA и др. Следует отметить, что применение аппроксимирующих функций формы непосредственно к компонентам вектора перемещения внутренней точки конечного элемента корректно только при выполнении расчетов в декартовой системе координат. При выполнении расчетов на основе МКЭ в криволинейной системе координат такой прием аппроксимации приводит к общеизвестной проблеме учета смещения конечного элемента как твердого тела [16; 23]. В настоящей работе для определения НДС нелинейно деформируемой оболочки вращения без использования гипотезы прямой нормали разработан на шаге нагружения шестигранный конечный элемент в формулировке метода перемещений с узловыми неизвестными в виде приращений перемещений и их производных. Для аппроксимации приращений перемещений через узловые значения разработана векторная аппроксимация класса C(1), использование которой после координатного преобразования дало возможность получить аппроксимирующие функции, приводящие к решению проблемы учета смещения конечного элемента как абсолютно твердого тела Методы Векторные параметры оболочки вращения Положение произвольной точки M0 срединной поверхности оболочки вращения в декартовой системе координат 0х, у, z с ортами определяется радиус-вектором (1) где r - радиус вращения точки M0. Векторы базиса произвольной точки M0 определяются выражениями Радиус-вектор произвольной точки оболочки имеет вид (2) Векторы базиса произвольной точки определяются дифференцированием (2) и представляются матричными выражениями (3) где Дифференцированием (2) с использованием (3) можно производные базисных векторов точки представить компонентами в базисе этой же точки: (4) где Перемещения и деформации В условиях шагового нагружения рассматриваются три положения произвольной точки оболочки: исходное состояние, деформированное состояние (вектор перемещения ) и соседнее с деформированным состоянием (вектор перемещения ). Компоненты векторов перемещений и определяются в базисе точки : (5) После дифференцирования (5) с учетом (4) получаются соотношения (6) где - являются функциями компонент векторов перемещений , соответственно и определяются выражениями такого вида: (7) Определяя положение точки радиус-вектором (8) его дифференцированием можно найти векторы локального базиса точки : (9) Компоненты тензора деформаций после j шагов нагружения определяются разностью компонент метрических тензоров в точках и [3]: (10) При учете (6) и (9) деформации запишутся выражениями (11) Компоненты тензора приращений деформаций на (j+1)-м шаге нагружения определяются разностью ковариантных компонент метрических тензоров точек и : (12) где - компоненты линейного тензора приращений деформаций и нелинейного тензора приращений деформаций. При учете (6) и (7) линейные части тензора приращений деформаций можно представить выражениями (13) где Компоненты линейного тензора приращений деформаций можно определить матричной зависимостью (14) где - матрица операторов соотношений (13). Зависимости напряжений от деформации Полные напряжения и их приращения в актуальном базисе точки Mζ определяются соотношениями [3; 12] (15) где - контравариантные компоненты тензоров напряжений и их приращений; - ковариантные компоненты тензоров деформаций и их линейных приращений; - параметры Ламе; - первые инварианты тензоров деформаций и их линейных приращений. Соотношения (15) представляются в матричном виде: (16) где Матрица жесткости конечного элемента в форме шестигранника Узловые координаты шестигранного конечного элемента принимаются в виде матриц-строк: (17) где - координаты узлов конечного элемента в глобальной системе координат; - узловые точки шестигранника. Для выполнения численного интегрирования по объему шестигранного конечного элемента он отображается на куб, с локальными координатами, изменяющимися в пределах -11 [12]. Координаты внутренней точки шестигранного дискретного элемента определяются через узловые значения на основе трилинейных функций в системе координат куба: (18) Прямые и обратные производные координат определяются дифференцированием (18). Скалярная аппроксимация перемещений. Принятые узловые неизвестные шестигранника в локальной и глобальной системах координат записываются матрицами-строками: (19) Матрицы-столбцы (19) связаны между собой матричной зависимостью (20) где матрица содержит производные глобальных координат узловых точек в локальной системе Обычно в МКЭ компоненты вектора аппроксимируются выражениями (21) где - аппроксимирующая матрица с элементами в виде полиномов Эрмита третьей степени; (22) После дифференцирования (21) находятся производные компонент вектора перемещения по глобальным координатам (23) В аппроксимирующих выражениях (21) и (23) отсутствуют параметры используемой для расчета криволинейной системы координат, поэтому каждая определяется через узловые значения только этой же компоненты. Указанные недостатки аппроксимирующих выражений (21) и (23) приводят в МКЭ к отсутствию возможности учета элемента как твердого тела. С учетом (21) и (23) соотношения (14) для приращений деформаций представляются в матричном виде (24) где - строка узловых неизвестных объемного конечного элемента; Векторная аппроксимация перемещений. При выполнении расчетов в криволинейных системах координат для получения аппроксимирующих выражений искомых величин вводят в рассмотрение векторные узловые неизвестные конечного элемента в локальной и глобальной системах координат выражениями (25) Между матрицами-столбцами (25) выполняется соотношение (26) Традиционная теория аппроксимаций позволяет определить вектор перемещения внутренней точки конечного элемента через узловые векторы (25) выражением (27) Дифференцированием (27) по криволинейным координатам определяются производные векторов перемещений: (28) На основании (5), (6) столбец узловых неизвестных можно представить матричным выражением (29) где (30) - векторы базисных узловых точек При координатном преобразовании выражения (29) базисные векторы узловых точек конечного элемента в матрице должны быть выражены через базисные векторы внутренней точки конечного элемента соотношением (31) Используя (31) для замены элементов матрицы в (29) и подставляя (29) в (27) и (28), можно сформировать матричные выражения (32) Приравнивая правые части (5), (6) и (32), можно получить аппроксимирующие выражения (33) где матрицы и являются строками матриц и соответственно. Компонентами матрицы в аппроксимациях (33) учитываются параметры используемой в расчете криволинейной системы координат. Отдельная компонента вектора перемещения зависит от всех компонент векторов перемещений узловых точек их производных, что более адекватно соответствует геометрическому смыслу в криволинейной системе координат. С учетом (33) компоненты линейного тензора приращений деформаций запишутся матричным выражением (34) Матрица жесткости конечного элемента на шаге нагружения. Используется функционал, основанный на равенстве работ внешних и внутренних сил на шаге нагружения. (35) где q и Δq - суммарная и шаговая нагрузка на поверхности s элемента; V - объем элемента. После подстановки аппроксимирующих выражений и минимизации функционала по узловым неизвестным получается выражение (36) где [K] - матрица жесткости элемента; [Kн] - матрица за счет нелинейной части приращений деформаций; - вектор узловых нагрузок на шаге нагружения; - невязка на шаге нагружения. Результаты и обсуждение Соотношениями (36) представляются два варианта матрицы жесткости шестигранного конечного элемента в криволинейной системе координат. Первый вариант основан на аппроксимации перемещений, корректной только в декартовой системе координат, когда каждая искомая величина аппроксимируется через узловые значения этой же величины. Второй вариант основан на математической модели аппроксимации векторных величин, в которой после координатного преобразования получаются аппроксимирующие выражения искомых величин, включающие параметры используемой криволинейной системы координат. Пример расчета 1. Определено напряженно-деформированное состояние цилиндрической оболочки диаметром d и длиной L, находящейся под действием сосредоточенной силы Р (рис. 1). Геометрические параметры оболочки приняты следующими: L = 0,8 м; d = 0,504 м; t = 0,0254 м; E = 6,9×104 МПа; n = 0,28; Р = 6,9 н. t Рис. 1. Цилиндрическая оболочка на пружинных опорах Figure 1. Cylindrical shell on spring supports В качестве опорных устройств в расчете приняты пружины переменной жесткости. В первом варианте расчета считалось, что пружины абсолютно жесткие. Во втором варианте предполагалось, что пружины позволяют смещение всей конструкции на величину Δ. Расчеты выполнялись с использованием разработанного шестигранного конечного элемента на основе двух вариантов аппроксимации искомых величин. В первом варианте использовалась традиционная аппроксимация компонент вектора перемещения (21) и их производных (23). Во втором варианте расчета применялась разработанная авторами аппроксимация перемещений и их производных (33). Результаты вариативных расчетов при абсолютной жесткости пружинных опор представлены в таблице в зависимости от густоты сетки дискретизации четвертой части оболочки. Результаты вариативных расчетов при абсолютной жесткости пружинных опор Густота сетки Прогиб w, м Вариант расчета 1 Вариант расчета 2 3×3×1 2,15·10-3 2,0·10-3 4×4×1 2,18·10-3 2,1·10-3 5×5×1 2,19·10-3 2,1·10-3 6×6×1 2,2·10-3 2,2·10-3 Results of variable calculations with absolute stiffness of spring supports Mesh density Deflection w, m Calculation option I Calculation option II 3×3×1 2.15·10-3 2.0·10-3 4×4×1 2.18·10-3 2.1·10-3 5×5×1 2.19·10-3 2.1·10-3 6×6×1 2.2·10-3 2.2·10-3 Табличные результаты свидетельствуют, что значения перемещения точки приложения силы Р совпадают в случае использования каждого из вариантов аппроксимации искомых величин шестигранника. На рис. 2 представлен график изменения окружных нормальных напряжений внутренних волокон точки Р в зависимости от жесткости пружинных опор, позволяющих смещение конструкции на величину Δ. Как видно из графика, смещение цилиндрической оболочки на 0,1 м как твердого тела не приводит к изменению окружных нормальных напряжений (линия II). Использование аппроксимации искомых величин МКЭ в варианте I приводит к значительным изменениям окружного нормального напряжения (линия I). Значения нормальных напряжений по варианту II остались неизменными и при Δ = 1 м, что свидетельствует о решении проблемы учета смещения конечного элемента как твердого тела на основе разработанного варианта II аппроксимации компонент вектора перемещения через узловые значения компонент. σк IIIΔ, м Рис. 2. График изменения окружного нормального напряжения Figure 2. Graph of changes in the circumferential normal stress Пример расчета 2. Рассматривалось деформированное состояние цилиндрической панели, защемленной по концам [22] (рис. 3). В качестве исходных данных приняты следующие величины: толщина поперечного сечения арки - t = 0,00476 м; ширина поперечного сечения арки - b = 0,00254 м; внутренний радиус R = 3,381 м; сектор круговой арки α = 0,256 рад.; модуль упругости материала арки Е = 7·104 МПа; коэффициент Пуассона ν = 0,2. Нагружение цилиндрической панели осуществлялась сосредоточенной силой Р, прикладываемой в ее вершине. На графике (рис. 4) приведены значения перемещений точки приложения груза в зависимости от величины силы, где использованы обозначения: Р, кН - значение сосредоточенной силы; w1, м - перемещения точки, найденные при использовании скалярной аппроксимации перемещений в дискретном элементе; w2, м - численные значения перемещений, найденные на основе аппроксимации перемещений в векторной формулировке; w3, м - перемещения, приведенные в [22] на основе итерационной процедуры. Анализ приведенных на графике численных результатов показывает, что значения перемещений w1, полученные на основе использования скалярной аппроксимации, оказались отличающимися от двух других в пределах 8 % при нагрузках меньших 0,1 кН. При большем значении нагрузки в случае скалярного варианта аппроксимации перемещений происходит сбой вычислительного процесса. Результаты перемещений, полученные на основе аппроксимации искомых величин в векторной формулировке, находятся в хорошем соответствии с результатами [22], что свидетельствует о корректности разработанного алгоритма векторной аппроксимации перемещений. w0 Рис. 3. Расчетная схема цилиндрической панели, защемленной на концах Figure 3. Design diagram of a cylindrical panel pinched at the ends w3w1w2w, мP, кН Рис. 4. График значений перемещений w точки приложения сосредоточенной силы P Figure 4. Graph of displacement values w of the point of application of concentrated force P Заключение Векторная аппроксимация искомых величин МКЭ при расчетах в криволинейных системах координат является более корректной, так как позволяет учитывать в аппроксимирующих выражениях параметры используемой криволинейной системы координат. При использовании скалярного варианта аппроксимации искомых величин МКЭ тип криволинейной системы координат во внимание не принимается, что приводит к некорректности аппроксимирующих соотношений.

Об авторах

Наталья Анатольевна Гуреева

Финансовый университет при Правительстве Российской Федерации

Email: nagureeve@fa.ru
ORCID iD: 0000-0003-3496-2008

доктор физико-математических наук, профессор, доцент департамента математики

Российская Федерация, 125993, Москва, Ленинградский пр-кт, д. 49

Румия Зайдуллаевна Киселева

Волгоградский государственный аграрный университет

Email: rumia1970@yandex.ru
ORCID iD: 0000-0002-3047-5256

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

Российская Федерация, 400002, Волгоград, Университетский пр-кт, д. 26

Анатолий Петрович Киселев

Волгоградский государственный аграрный университет

Email: apkiselev1969@yandex.ru
ORCID iD: 0000-0002-7138-2056

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

Российская Федерация, 400002, Волгоград, Университетский пр-кт, д. 26

Анатолий Петрович Николаев

Волгоградский государственный аграрный университет

Email: anpetr40@yandex.ru
ORCID iD: 0000-0002-7098-5998

доктор технических наук, профессор кафедры механики, инженерно-технологический факультет

Российская Федерация, 400002, Волгоград, Университетский пр-кт, д. 26

Юрий Васильевич Клочков

Волгоградский государственный аграрный университет

Автор, ответственный за переписку.
Email: klotchkov@bk.ru
ORCID iD: 0000-0002-1027-1811

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

Российская Федерация, 400002, Волгоград, Университетский пр-кт, д. 26

Список литературы