# Abstract

Introduction. With the development of numerical methods and computational complexes, it is quite easy to evaluate the stress state of thin-walled structures in the form of rotation bodies. However, when solving such problems by the finite element method, it is necessary to choose such finite element grid to "grasp" all possible singularities of the stressed state. To correctly take them into account, you must reduce the size of the finite elements. Reducing the size of the elements leads to an increase in the required computing power.

Formulation of the problem. When solving applied problems, even with a sufficiently coarse grid, the number of elements can exceed hundreds of thousands. When solving problems for real constructions in a three-dimensional setting, the amount of computation can be quite large and not every supercomputer can even handle such a solution.

Objective. The purpose of this paper is to use the well-known approach used in shell theory, which allows us to reduce the three-dimensional problem to the solution of a one-dimensional problem, which substantially reduces the requirements for computing power.

Method (methodology). The problem of determining the stress state of shell structures in the form of bodies of revolution is considered. The approach is based on the integration of the equations of the theory of shells and the expansion of functions into Fourier series for separation of variables. The expansion into a discrete Fourier series in cosines and sines is used in this paper, which describes arbitrary asymmetric mechanical loads.

Results. A thin-walled cylindrical structure hinged at the ends is considered. The structure is loaded in three places by a distributed force acting normal to the surface of the shell. After integrating the system of equations for the shell, the found stress-strain state of the shell is determined by the stress components on the outer and inner surfaces of the shell and the displacement components. The paper compares the calculation results with the proposed methodology and the finite element method.

The conclusion. It is shown that the use of methods of shell theory, and the proposed expansion of resolving functions and loads in a Fourier series, allows solving problems using small computing resources. At the same time, the necessary accuracy of calculation for all components of the stress-strain state of the structure is ensured.

# Full Text

## Введение

В настоящее время при определении напряженного состояния различных конструкционных элементов широко используются вычислительные комплексы, основанные на методе конечных элементов. Произвести оценку напряженного состояния тонкостенных конструкций, в виде тел вращения используя такие комплексы как ANSYS, COSMOS [1, 2] и др. не представляет труда даже в геометрически и физически нелинейной постановке. Напряженное состояние может быть найдено в одномерной, двумерной или трехмерной постановке в зависимости от типа и цели задачи и располагаемой вычислительной мощности. Однако существует класс задач, где заложенная в данных вычислительных комплексах универсальность, теряет свою эффективность. Это происходит при решении задач с подвижными границами, т.е. контактных задач для конструкций, взаимодействующих с жесткими или упругими телами. При решении задач методом конечных элементов для тонкостенных конструкций, взаимодействующих с другими телами, необходимо выбрать такую сетку конечных элементов, чтобы «схватить» все возможные сингулярности напряженного состояния, но при этом использовать минимально возможную размерность задачи. Этого не всегда легко добиться, поскольку, например, при контакте оболочечных конструкций возникают локализованные нагрузки на границе области контакта [3, 4, 5]. Чтобы корректно их учитывать, необходимо уменьшать размер конечных элементов. Уменьшение размеров элементов приводит к увеличению требуемых вычислительных мощностей. При решении реальных прикладных задач даже при достаточно грубой сетке количество элементов может превышать сотни тысяч. Например, при определении не осесимметричного напряженного состояния (двумерная задача) оболочки вагона-цистерны потребовалось около 100 тыс. оболочечных конечных элементов типа shell [6]. При решении такой задачи в одномерной постановке (осесимметричное напряженное состояние) необходимо около 1 тыс. элементов. При решении данной задачи в трехмерной постановке (в рамках теории упругости) необходимо не менее  объемных конечных элементов типа solid для случая если толщина оболочки будет аппроксимирована одним конечным элементом (при необходимом соблюдении соотношений между размерами конечного элемента). Если для увеличения точности расчета толщина оболочки будет аппроксимирована тремя элементами, то число элементов возрастет до конечных элементов. Если толщина оболочки будет аппроксимирована шестью элементами необходимо не менее  элементов. Если для первого расчетного случая для получения решения необходимо около суток, то для третьего случая речь пойдет уже о месяцах, к тому же объем вычислений будет достаточно велик и далеко не всякий компьютер сможет такое решение обработать. Эти оценки приведены для решения задачи для оболочки вагона-цистерны на суперкомпьютере «Уран» ИММ УрО РАН. Следовательно, решение такой задачи в трехмерной постановке становится затруднительным даже при использовании суперкомпьютеров. Поэтому вопрос об уменьшении размерности задачи для решения прикладных задач для тонкостенных конструкций является весьма актуальным.

## Постановка задачи

Одним из эффективных методов понижения размерности задачи для тонкостенных конструкций в виде тел вращения является использования методов теории оболочек и применение рядов Фурье. Известно, что в методе конечных элементов для тел вращения при аппроксимации внешних нагрузок также используют разложение в ряд Фурье по окружной координате [7, 8]. Однако при этом размерность задачи не понижается.

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

Если для описания тонкостенной конструкции использовать оболочечные модели, основанные на различных гипотезах, то напряженное состояние в ортогональной криволинейной системе координат, описывается системой дифференциальных уравнений в частных производных [10, 11, 12]

(1)

Здесь n –порядок производных, зависящий от принятых допущений; ${A}_{m}$ - матрицы, элементы которых определяются через геометрические и жесткостные параметры оболочки; $\overline{f}$ - вектор, компоненты которого зависят от приложенных к оболочке поверхностных нагрузок и интегральных характеристик температурного поля. Для классической теории оболочек искомый вектор $\overline{N}$ имеет восьмой порядок, а для уточненных моделей – десятый (модель с учетом гипотезы Тимошенко), двенадцатый (модель с учетом сдвига и обжатия) и др.

Система (1) дополняется граничными условиями на контуре $\alpha =const$. Поскольку рассматриваются оболочки вращения замкнутые в направлении $O\beta$, то граничные условия в этом направлении заменяются условиями периодичности. Краевая задача системы (1) за счет периодичности позволяет для всех искомых функций представить решения в виде рядов Фурье по координате $\beta$.

## Материалы и методы исследования

В данной работе будем использовать классическую теорию оболочек, основанную на гипотезах Кирхгофа-Лява [10]. Координатную поверхность оболочки вращения отнесем к криволинейной ортогональной системе $s$,$\theta$. На рис. 1 показаны: $s$ - длина дуги меридиана, $\theta$ - центральный угол в круге, перпендикулярном оси вращения z. Следовательно, задача определения напряженного состояния оболочки сводится к системе [13]

(2)

$\overline{Y}=${},   (3)

где $\overline{Y}$ - вектор разрешающих функций;  - радиальное и осевое усилия;  - радиальное и осевое перемещения; $\stackrel{^}{S}$ - сдвигающее усилие; ${M}_{s}$ - меридиональный изгибающий момент; $v$ - окружное перемещение; ${\vartheta }_{s}$ - угол поворота нормали; Е - единичная матрица, $\Omega$ - область контакта.

В отличие от уравнения (1) в правую часть уравнения (2) добавлено еще одно слагаемое, которое учитывает действие контактной нагрузки ${q}_{c}$ (по радиальному, осевому и окружному направлению), если оболочка контактирует с основанием.

Внешнюю нагрузку, действующею на оболочку вращения можно представить в виде компонент распределенной нагрузки действующих по касательной к образующей ${q}_{s}$ и направляющей ${q}_{\theta }$ и по нормали ${q}_{\gamma }$ к поверхности оболочки. Поскольку обычно для оболочки вращения легко можно подобрать поперечную ось, относительно которой внешние нагрузки будут симметричны или антисимметричны, то их компоненты можно представить в виде разложения [7, 10].

${q}_{s}\left(s,\theta \right)=\sum _{k=0}^{\infty }{q}_{s}\left(s\right)\mathrm{cos}k\theta$,

${q}_{\theta }\left(s,\theta \right)=\sum _{k=1}^{\infty }{q}_{\theta }\left(s\right)\mathrm{sin}k\theta$, (4)

${q}_{\gamma }\left(s,\theta \right)=\sum _{k=0}^{\infty }{q}_{\gamma }\left(s\right)\mathrm{cos}k\theta$,

где $k$ - номер гармоники.

В силу периодичности компонент поверхностной нагрузки все функции вектора $\overline{Y}$ можно разложить в ряд Фурье по окружной координате $\theta$. Для симметричных компонент вектора  имеем

${\Phi }^{sim}\left(s,\theta \right)=\sum _{k=0}^{\infty }{\Phi }_{i}^{sim}\left(s\right)\mathrm{cos}k\theta$ (5)

а, для антисимметричных компонент ,

${\Phi }^{ans}\left(s,\theta \right)=\sum _{k=1}^{\infty }{\Phi }_{i}^{ans}\left(s\right)\mathrm{sin}k\theta$ (6)

Индексом sim отмечены симметричные а, индексом ans – антисимметричные компоненты вектора.

Учитывая (4) - (6) система в частных производных (2) сводится к ряду систем обыкновенных дифференциальных уравнений в нормальной форме восьмого порядка [10]

(7)

с граничными условиями

${B}_{1}\overline{Y}\left({s}_{0}\right)={\overline{b}}_{1}$

Здесь  - заданные матрицы, - заданные векторы. Вектор нагрузок $\overline{{g}_{k}}$ состоит из компонентов внешней и контактной нагрузки. Компоненты контактной нагрузки должны быть заранее найдены из решения контактной задачи. Методы решения одномерных и двумерных контактных задач и нахождение вектора $\overline{{g}_{k}}$ подробно рассмотрены в работе [13].

Для решения краевой задачи системы (7) будем использовать метод дискретной ортогонализации [14]. Поскольку система (7) содержит амплитудные значения разрешающих функций, то все нагрузки и функции (4) - (6) будем строить на дискретном множестве точек, т.е. применять дискретные ряды Фурье [15].

Таким образом, для интегрирования уравнения (7) учитывая (4) значение нагрузки должно быть представлено в виде разложения в дискретный ряд Фурье [15, 16].

(8)

Поскольку используем подход, основанный на применении дискретных рядов Фурье для функций, то для этого на поверхность оболочки нанесем криволинейную сетку с равным шагом по меридиану и окружности. Таким образом, получим множество виртуальных оболочечных элементов. На рис. 1 показан виртуальный элемент с размерами  и ${a}_{\theta }\text{\hspace{0.17em}}$ по меридиану и окружности, нагруженный распределенной нагрузкой ${q}_{\gamma }$. Таким образом, будем считать, что на любом виртуальном элементе известно значение компонент внешней и контактной нагрузки.

Рис. 1. Цилиндрическая оболочка, нагруженная внешней нагрузкой q распределенной на виртуальном элементе

При решении задач с использованием теории оболочек [10] обычно рассматриваются нагрузки симметричные относительно поперечной оси. Следовательно, при вертикальной оси симметрии X, показанной на рис.1, можно использовать разложение (4).

При решении некоторых задач для тел вращения могут встречаться нагрузки не симметричные относительно поперечной оси. Примером таких нагрузок могут быть ветровые нагрузки, силы инерции при землетрясениях, когда земля получает ускорение и другие виды нагрузок. Следовательно, необходимо аппроксимировать нагрузки, распределенные произвольным образом, которые можно представлять в виде [8].

${q}_{s}\left(s,\theta \right)=\sum _{k=0}^{\infty }{q}_{s}^{sim}\left(s\right)\mathrm{cos}k\theta +\sum _{k=1}^{\infty }{q}_{s}^{ans}\left(s\right)\mathrm{sin}k\theta$,

${q}_{\theta }\left(s,\theta \right)=-\sum _{k=1}^{\infty }{q}_{\theta }^{sim}\left(s\right)\mathrm{sin}k\theta +\sum _{k=0}^{\infty }{q}_{\theta }^{ans}\left(s\right)\mathrm{cos}k\theta$,                                                             (9)

${q}_{\gamma }\left(s,\theta \right)=\sum _{k=0}^{\infty }{q}_{\gamma }^{sim}\left(s\right)\mathrm{cos}k\theta +\sum _{k=1}^{\infty }{q}_{\gamma }^{ans}\left(s\right)\mathrm{sin}k\theta$.

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

${\Phi }^{ans}\left(s,\theta \right)=\sum _{k=1}^{\infty }{\Phi }_{i}^{ans}\left(s\right)\mathrm{sin}k\theta$, (10)

и, для антисимметричных компонент $\stackrel{^}{S}$$v$

${\Phi }^{sim}\left(s,\theta \right)=\sum _{k=0}^{\infty }{\Phi }_{i}^{sim}\left(s\right)\mathrm{cos}k\theta$. (11)

Следовательно, напряженно-деформированное состояние оболочки будет определяться суммой компонент вектора разрешающих функций (3)

$\overline{Y}={\overline{Y}}_{1}+{\overline{Y}}_{2}$, (12)

где ${\overline{Y}}_{1}$ - вектор разрешающих функций при разложении (5), (6), а ${\overline{Y}}_{2}$ при разложении (10), (11).

Таким образом, поскольку функция нагрузки известна на множестве виртуальных элементов, периодическая с периодом N, то ее необходимо разложить в ряд Фурье по косинусам и синусам. Среди разнообразных разложений в ряд Фурье дискретных функции будем использовать коэффициенты, дающие хорошую точность аппроксимации [17] .

(13)

,

где ${q}_{i}$ - значение нагрузки на каждом виртуальном элементе.

Для демонстрации рассмотренного выше подхода рассмотрим тонкостенную цилиндрическую конструкцию, шарнирно закрепленную по торцам. Конструкция нагружена в трех местах распределенной силой $q$, действующей по нормали к поверхности оболочки. На рис. 1. показан один нагруженный виртуальный элемент.

При расчете принималось - длина оболочки L=2.40 м, внешний радиус ${R}_{}$= 1,50 м, толщина h = 0,006 м, модуль упругости $E=2,1×{10}^{5}$МПа, коэффициент Пуассона n = 0,3. Граничные условия приняты

${N}_{r}={N}_{z}={M}_{s}=\stackrel{^}{S}=0,$ при $s=0$,

${N}_{r}={N}_{z}={M}_{s}=\stackrel{^}{S}=0$, при $s=L$.

Окружность разбивалась углом равным q =2°, который стягивает виртуальный элемент длиной ${a}_{\theta }\text{\hspace{0.17em}}$, т.е. длина окружности оболочки аппроксимирована 180 виртуальными элементами. Длина оболочки аппроксимирована 10 виртуальными элементами с длиной =0.24 м.

Рассмотрен случай, когда на оболочку действуют три силы ${P}_{1}$=1000 Н, ${P}_{2}$=1000 Н и ${P}_{3}$=3000 Н, равномерно распределенных на трех виртуальных элементах. Координаты центра нагруженных элементов $s$=L/2 и $\theta$=$\pi$/3, $s$=L/2 и $\theta$=$\pi$/2, $s$=L/2 и $\theta$=2$\pi$/3 соответственно. При действии такой нагрузки не представляется возможным подобрать поперечную ось оболочки, чтобы использовать стандартную процедуру разложения в ряд (4) при действии симметричных нагрузок. Поэтому для решения этой задачи необходимо использовать суперпозицию двух решений (12), рассмотренную выше. При этом численным экспериментом установлено, что в разложении в ряд Фурье (8) необходимо удерживать $k=N/2$ гармоник.

После интегрирования системы (7) найденное напряженно-деформированное состояние оболочки определяется компонентами напряжений ( ${\sigma }_{s}$- меридиональное, ${\sigma }_{\theta }$ -окружное, ${\tau }_{s\theta }$ -касательное напряжение) на внешней и внутренней поверхностях оболочки и компонентами перемещений срединной линии [10].

Точность решения задачи зависит от используемой модели оболочки (1), от точности интегрирования системы уравнений (7) и от точности аппроксимирующей нагрузки (13). Эта точность может быть определена сопоставлением с решением аналогичной задачи другим численным методом или с экспериментом. В настоящее время, как отмечено выше, наиболее универсальным численным методом исследования тонкостенных конструкции является метод конечных элементов. Поэтому результаты различных компонентов напряжеpнно-деформированного состояния для исследуемой оболочки сравнивались с результатами, полученными с помощью вычислительного комплекса ANSYS [1], основанного на методе конечных элементов.

## Результаты

На рис.2 сплошной линией показано изменение по окружности радиального перемещения оболочки ${u}_{r}$ в зоне приложения сил. На рис.2 штриховой линией показано изменение радиального перемещения ${u}_{r}$, полученное при использовании программы ANSYS. Максимальное отклонение результатов расчета, вычисленное как $\delta =|\Delta {u}_{r}|×{{u}_{rcр}}^{-1}×100%$ (где ${u}_{rcр}$ - среднее значение перемещения) не превышает 6.7%. На рисунке видны три характерных локальных прогиба оболочки в местах приложения нагрузки при $\theta$=$\pi$/3, $\theta$=$\pi$/2, $\theta$=2$\pi$/3. При нагрузке ${P}_{3}$ прогиб, естественно, увеличенный.

Рис.2. Изменение радиального перемещения оболочки

Рис. 3. Изменение меридионального напряжения ${{\sigma }_{s}}^{-}$ на внутренней поверхности оболочки

На рис.3 сплошной линией показано изменение меридионального напряжения по окружности на внутренней поверхности оболочки ${{\sigma }_{s}}^{-}$ в зоне приложения сил. На рис.3 штриховой линией показано изменение меридионального напряжения, полученное при использовании программы ANSYS. Максимальное отклонение результатов расчета, вычисленное как $\delta =|\Delta {\sigma }_{s}^{-}|×{{{\sigma }_{s}^{-}}_{cр}}^{-1}×100%$ (где ${{\sigma }_{s}^{-}}_{cр}$- среднее значение напряжения) не превышает 5.5%.

На рис.4 сплошной линией показано изменение окружного напряжения на внутренней поверхности оболочки ${{\sigma }_{\theta }}^{-}$ по окружности в зоне приложения сил. На рис.4 штриховой линией показано изменение окружного напряжения на внутренней поверхности оболочки, полученного при использовании программы ANSYS. Максимальное отклонение $\delta =|\Delta {\sigma }_{\theta }^{-}|×{{{\sigma }_{\theta }^{-}}_{cр}}^{-1}×100%$ (где ${{\sigma }_{\theta }^{-}}_{cр}$ - среднее значение напряжения) не превышает 0.6%.

Рис. 4. Изменение окружного напряжения на внутренней поверхности оболочки

Из рисунков 2 - 4 видно, что характер напряженно-деформированного состояния исследуемой оболочки вычисленный, с использованием двух различных численных методов полностью совпадает.

## Обсуждение результатов

При расчете методом конечных элементов окружность аппроксимирована 180 четырех узловыми элементами shell-181. Поскольку точки выдачи напряженного состояния в методе конечных элементов предусмотрена по узлам, а в предложенной методе в середине элемента, то происходит сдвижка в представлении информации, что приводит к расхождению результатов в некоторых точках до 6.7%. Следует так же отметить, что данное расхождение находится в пределах допустимой погрешности (7-8%), характерной для метода конечных элементов [1].

## Заключение

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

### Igor G. Emel’yanov

Federal State Budgetary Scientific Institution “Institute of Engineering Science, Ural Branch of the Russian Academy of Sciences”; Ural State University of Railway Transport

Author for correspondence.
Email: emelyanov@imach.uran.ru

Russian Federation

chief research officer of Institute of Engineering Science; professor

### Alexey V. Kuznetsov

Federal State Budgetary Scientific Institution “Institute of Engineering Science, Ural Branch of the Russian Academy of Sciences”

Email: Alekseikuz031082@rambler.ru

Russian Federation

Research Fellow

# References

1. Басов К. А. ANSYS: Справочник пользователя / К. А. Басов. - М.: ДМК Пресс, 2005. - 640 с.
2. Алямовский А. А. SolidWorks / COSMOSWorks. Инженерный анализ методом конечных элементов. - М.: ДМК Пресс, 2004. - 432 с.
3. Григолюк Э. И. Контактные задачи теории пластин и оболочек / Э. И. Григолюк, В. М. Толкачев. - М.: Машиностроение, 1980. - 411 с.
4. Джонсон К. Механика контактного взаимодействия / К. Джонсон. - М.: Мир, 1989. - 510 с.
5. Моссаковский В. И. Контактное взаимодействие элементов оболочечных конструкций / В. И. Моссаковский, В. С. Гудрамович, Е. М. Макеев. - Киев: Наукова думка, 1988. - 288 с.
6. Емельянов И. Г. Определение напряженного состояния и ресурса оболочечной конструкции / И. Г. Емельянов, В. И. Миронов, А. В. Кузнецов // Проблемы машиностроения и надежности машин. - 2007. - №5. - С. 57-65.
7. Зенкевич О. К. Метод конечных элементов в теории сооружений и механике сплошных сред / О. К. Зенкевич, И. Чанг. - М.: Недра, 1974. - 240 с.
8. Галлагер Р. Метод конечных элементов. Основы / Р. Галлагер. - М.: Мир, 1984. - 428с.
9. Емельянов И. Г. Применение виртуальных элементов при определении напряженного состояния оболочек вращения / И. Г. Емельянов, А. В. Кузнецов // Вычислительная механика сплошных сред. - 2014. - Т. 7. - № 3. - С. 245-252.
10. Григоренко Я. М. Методы расчета оболочек: Теория оболочек переменной жесткости / Я. М. Григоренко, А. Т. Василенко. - Киев: Наукова думка, 1981. - Т. 4. - 544 с.
11. Григоренко Я. М. Численно-аналитическое решение задач механики оболочек на основе различных моделей / Я. М. Григоренко, Г. Г. Влайков, А. Я. Григоренко. - Киев: Академпериодика, 2006. - 472 с.
12. Василенко А. Т. Расчет параметров напряженного состояния конструктивных элементов из композиционных материалов на основе оболочечных моделей / А. Т. Василенко, Г. П. Голуб, Я. М. Григоренко // Расчеты на прочность. Вып. №30. - М.: Машиностроение, 1989. - С. 87-96.
13. Емельянов И. Г. Контактные задачи теории оболочек / И. Г. Емельянов. - Екатеринбург: УрО РАН , 2009. - 185с.
14. Годунов С. К. О численном решении краевых задач для систем линейных обыкновенных дифференциальных уравнений / С. К. Годунов // Успехи математических наук. - 1961. - Вып. 16. - № 3. - С. 171-174.
15. Григоренко Я. М. Решение краевых задач о напряженном состоянии упругих тел сложной геометрии и структуры с применением дискретных рядов Фурье / Я. М. Григоренко // Прикладная механика. - 2009. - Вып. 45. - №5. - С. 3-51.
16. Корн Г. Справочник по математике для научных работников и инженеров / Г. Корн, Т. Корн. - М.: Наука, 1968. - 720 с.
17. Ба-хуссейн А. А. Дискретное преобразование Фурье. - URL: http://ilab.xmedtest.net/?q=node/3740 (01.02.2015).

# Statistics

#### Views

Abstract - 412

PDF (Russian) - 306

#### PlumX

Copyright (c) 2017 Emel’yanov I.G., Kuznetsov A.V.