Determination of the stressed state of thin-construction structures using the methods of the theory of shells

Full Text

Введение

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

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

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

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

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

$\frac{\partial \overline{N}}{\partial \alpha }=\sum _{m=0}^n{A}_m\left(\alpha \right)\frac{{\partial }^m\overline{N}}{\partial {\beta }^m}+\overline{f}\left(\alpha ,\beta \right)\text{ .}$ (1)

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

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

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

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

$\frac{\partial \overline{Y}}{\partial s}=\sum _{m=0}^4{A}_m\left(s\right)\frac{{\partial }^m\overline{Y}}{\partial {\theta }^m}+\overline{f}\left(s,\theta \right)+\lambda E{\overline{q}}_c\left(s,\theta \right)\text{ ,}$ (2) 

$\overline{Y}=${$N_r, N_z, \hat{S}, M_s, u_r, u_z, v, {\vartheta }_s$}, ${}_{\text{c}}\text{=}{\left\{{\text{q}}_{\text{r}}{\text{, q}}_{\text{z}},{\text{ q}}_{\theta },\text{ 0, 0, 0, 0, 0}\right\}}^T\text{ ,}$ $\lambda \left(\text{s,}\theta \in \Omega \right)=1,\lambda \left(\text{s,}\theta \notin \Omega \right)=0,$ (3)

где $\overline{Y}$ - вектор разрешающих функций; $N_r,\text{ }{N}_Z\text{ }$ - радиальное и осевое усилия; $u_r,u_z$ - радиальное и осевое перемещения; $\hat{S}$ - сдвигающее усилие; $M_s$ - меридиональный изгибающий момент; $v$ - окружное перемещение; ${\vartheta }_s$ - угол поворота нормали; Е - единичная матрица, $\Omega$ - область контакта.

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

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

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

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

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

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

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

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

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

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

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

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

$\begin{array}{l}\frac{d\overline{Y}}{ds}=A_k\left(s\right){\stackrel{-}{Y}}_k+{\overline{g}}_k\left(s\right),\left(k=0,1,2\dots \right),\\ A_{{\rm K}}\left(s\right)=||{\text{a}}_{\text{ij}}^{\left(k\right)}\left(s\right)||,\left(i,j=1,\dots 8\right),s_0\le \text{s}\le {\text{s}}_{\text{L}}\text{, }{\overline{g}}_k=\left\{g_{1k},\dots g_{8k}\right\},\end{array}$ (7)

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

 $B_1\overline{Y}\left(s_0\right)={\overline{b}}_1$

$B_2\overline{Y}\left(s_L\right)={\overline{b}}_2\text{ .}$                                                                                          

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

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

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

$q_f=\frac{a_0}{2}+\sum _{k=1}^N\left(a_k\cos k\theta +b_k\sin k\theta \right)$ (8)

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

 

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

 

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

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

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

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

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

Разрешающие функции должны при этом быть дополнены слагаемыми - для симметричных компонент $N_r, N_z, M_s, u_r, u_z, {\vartheta }_s$

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

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

${\Phi }^{sim}(s,\theta )=\sum _{k=0}^{\infty }{\Phi }_i^{sim}\left(s\right)\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] .

$a_k=\frac{1}{N}\sum _{i=\text{0}}^{\text{N-1}}{\text{q}}_{\text{i}}\cos \frac{2\pi ki}{N}\text{ , }$  (13)                      

$b_k=\frac{1}{N}\sum _{i=\text{0}}^{\text{N-1}}{\text{q}}_{\text{i}}\sin \frac{2\pi ki}{N}\text{ , k}=\text{0, }...\frac{\text{N}}{\text{2}}$ ,

где $q_i$ - значение нагрузки на каждом виртуальном элементе.

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

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

$N_r=N_z=M_s=\hat{S}=0,$ при $s=0$,

$N_r=N_z=M_s=\hat{S}=0$, при $s=L$.

Окружность разбивалась углом равным q =2°, который стягивает виртуальный элемент длиной $a_{\theta }$, т.е. длина окружности оболочки аппроксимирована 180 виртуальными элементами. Длина оболочки аппроксимирована 10 виртуальными элементами с длиной $a_s$=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 =\left|\Delta u_r\right|\times {u_{rcр}}^{-1}\times 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 =\left|\Delta {\sigma }_s^{-}\right|\times {{{\sigma }_s^{-}}_{cр}}^{-1}\times 100\%$ (где ${{\sigma }_s^{-}}_{cр}$- среднее значение напряжения) не превышает 5.5%.

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

 

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

 

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

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

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

Заключение

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

About the authors

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

chief research officer of Institute of Engineering Science; professor

Russian Federation

Alexey V. Kuznetsov

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

Email: Alekseikuz031082@rambler.ru

Research Fellow

Russian Federation

References

Supplementary files