Известия РАН. Механика твердого тела, 2022, № 3, стр. 142-152

Введение. Актуальной проблемой является решение задач механики деформирования элементов пластин и оболочек [1–4] в тех случаях, когда не удается методом Фурье разделить переменные в уравнениях с частными производными и упростить задачу решением обыкновенных дифференциальных уравнений [5, 6]. Проблема преодолевается решением уравнений с частными производными с контролируемой погрешностью. Закладываются основы построения конечных элементов оболочек, на которых математические модели их деформирования решаются с контролируемой погрешностью. Таким образом, становится возможным построения алгоритмов решения задач механики деформирования тонкостенных конструкций с определением мест концентрации напряжений и их максимальных значений, что необходимо при решении задач прочности, например, летательных аппаратов. В работе представлена методика решения краевых задач механики деформирования пологой оболочки двоякой кривизны. Методика основана на решении линейных обыкновенных дифференциальных уравнений удовлетворяющих произвольным начальным условиям. Вероятно, впервые такое решение было получено А.Н. Крыловым методом Коши [7]. Декомпозицией [8] уравнения с частными производными механики деформирования пологой оболочки приводятся к обыкновенным и алгебраическим уравнениям их связи. Решения однородных обыкновенных дифференциальных уравнений и частные решения для их правых частей получены на основе известных работ [9–12] в виде формул. Они опубликованы в работах [13–15]. Построен и реализован алгоритм решения на ЭВМ краевой задачи механики деформирования элемента пологой защемленной по краям оболочки при поверхностном давлении.

1. Реализация. Рассматривается пологая оболочка двоякой кривизны, прямоугольная в плане, контур которой образован меридианами и параллелями, параллельными линиям главных кривизн (рис. 1). В качестве математической модели механики деформирования оболочки используются общие разрешающие уравнения оболочки двоякой кривизны в перемещениях:

где: u, ${v}$, w – компоненты перемещений, ${{q}_{1}}$, ${{q}_{2}}$, ${{q}_{3}}$ – соответствующие перемещениям компоненты внешних поверхностных сил, r₁, r₂ – главные радиусы кривизны оболочки, h и μ – толщина оболочки и коэффициент Пуассона, B и D – мембранная и изгибная жесткости.

Методом декомпозиции [8] система уравнений в частных производных (1.1) приводится к следующей системе дифференциальных уравнений:

и уравнений их связи:

В системе уравнений (1.2) $f_{1}^{{\left( 1 \right)}}(x,y) \ldots f_{3}^{{\left( 3 \right)}}(x,y)$ – неизвестные функции, введенные при декомпозиции. Эти функции подлежат определению путем удовлетворения уравнений связи (1.3). Проанализируем уравнения системы (1.2).

Первое, четвертое и седьмое и уравнение системы (1.2) – уравнения в частных производных, содержащие производные только по x. Эти уравнения будем рассматривать как обыкновенные дифференциальные уравнения (ОДУ) с параметром y.

Второе, пятое и восьмое уравнение системы (1.2) – уравнения в частных производных, содержащее производные только по y. Аналогично, рассматриваем эти уравнения как ОДУ с параметром х.

Третье, шестое и девятое уравнение – уравнения со смешанными производными. Эти уравнения будем рассматривать как ОДУ с параметром, но не относительно искомой функции, а относительно одной из производных. Так на примере третьего уравнения $\frac{{1 + {{\mu }}}}{2}\frac{{{{\partial }^{2}}{v}_{1}^{{\left( 3 \right)}}}}{{\partial x\partial y}} = f_{1}^{{\left( 3 \right)}}(x,y)$ записано как $\frac{{\partial \left[ {\frac{{\partial {v}_{1}^{{\left( 3 \right)}}}}{{\partial y}}} \right]}}{{\partial x}} = \frac{2}{{1 + {{\mu }}}}f_{1}^{{\left( 3 \right)}}(x,y)$. Оно будет решено относительно $\frac{{\partial {v}_{1}^{{\left( 3 \right)}}}}{{\partial y}}$. В свою очередь ${v}_{1}^{{\left( 3 \right)}}$ определять из этого уравнения не будем, знания лишь величины $\frac{{\partial {v}_{1}^{{\left( 3 \right)}}}}{{\partial y}}$ достаточно для предлагаемого удовлетворения уравнениям связи.

Разобъем сторону оболочки вдоль оси Ox на m + 1 интервалов точками x₁₀, x₂₀, …, x_m0, а сторону оболочки вдоль оси Oy на n + 1 интервалов точками y₀₁, y₀₂, … , y_0n (рис. 1). Эти точки для значений параметров указанных выше ОДУ.

Далее нанесем сетку меридианов и параллелей через выбранные точки. Вдоль каждого меридиана и параллели теперь можно записать ОДУ с дифференцированием по $x$ или по $y$ соответственно. Меридианы и параллели пересекаются в точках. В точках пересечения меридианов и параллелей будем выполнять уравнения связи (1.3), в точках контура оболочки – краевые условия.

Для каждой точки (на контуре и на пересечении меридианов и параллелей) запишем по столбцу, содержащему перемещения и производные от них:

Комбинированием производных, входящих в столбцы (1.4) и (1.5), по известным формулам могут быть вычислены силовые факторы, характеризующие состояние сечения оболочки: нормальные, обобщенные сдвиговые силы, изгибающие моменты, обобщенные перерезывающие силы.

Неизвестные функции $f_{1}^{{\left( 1 \right)}}(x,y) \ldots f_{3}^{{\left( 3 \right)}}(x,y)$ аппроксимируем кусочно-постоянными. Для этого другими меридианами и параллелями разбиваем оболочку на области (рис. 2) таким образом, чтобы каждая из точек (например, 11, 12, 21, …, mn) попадала внутрь одной из областей. Принимаем значения неизвестных функций $f_{1}^{{\left( 1 \right)}}(x,y) \ldots f_{3}^{{\left( 3 \right)}}(x,y)$ постоянными и равными их значениям в точках 11, 12, 21, …, mn. Для каждой из точек 11, 12, 21, …, mn запишем столбец значений неизвестных функций:

Вдоль каждого меридиана и параллели запишем соответствующие дифференциальные уравнения системы (1.2) в матричной форме:

где: ${{a}_{{x27}}} = - \left( {\frac{1}{{{{r}_{1}}}} + \frac{{{\mu }}}{{{{r}_{2}}}}} \right)$, ${{a}_{{x92}}} = - \frac{{12}}{{{{h}^{2}}}}\left( {\frac{1}{{{{r}_{1}}}} + \frac{{{\mu }}}{{{{r}_{2}}}}} \right)$, ${{a}_{{x96}}}$ = $ - \frac{{12}}{{{{h}^{2}}}}\left( {\frac{1}{{r_{1}^{2}}} + \frac{1}{{r_{2}^{2}}} + \frac{{2{{\mu }}}}{{{{r}_{1}}{{r}_{2}}}}} \right)$, ${{c}_{{x44}}} = \frac{2}{{1 - \mu }}$, c_x53 = = $\frac{2}{{1 + \mu }}$

В краткой форме:

где: ${{a}_{{y57}}} = - \left( {\frac{1}{{{{r}_{1}}}} + \frac{{{\mu }}}{{{{r}_{2}}}}} \right)$, ${{a}_{{y95}}} = - \frac{{12}}{{{{h}^{2}}}}\left( {\frac{1}{{{{r}_{2}}}} + \frac{{{\mu }}}{{{{r}_{1}}}}} \right)$, ${{c}_{{y22}}} = \frac{2}{{1 - {{\mu }}}}$, ${{c}_{{y36}}} = \frac{2}{{1 + {{\mu }}}}$

Для каждой из точек 11, 12, 21, …, mn запишем в матричной форме уравнения связи:

Решения ОДУ (1.8), (1.10) определяются аналитически по формулам [9–11]. Вдоль меридиана y = y_0j между любыми соседними точками i, j и i, j – 1:

где: ${{\Delta }}x = {{x}_{{ij}}} - {{x}_{{i,j - 1}}}$ – расстояние по меридиану между соседними точками, ${\mathbf{K}}_{{i,j - 1}}^{{ij}}$ = = $\sum\nolimits_{k = 0}^\infty {\frac{{{{{({{{\mathbf{A}}}_{x}}\Delta x)}}^{k}}}}{{k!}}} $ – матрица функций Коши–Крылова (решение однородного ОДУ),

Вдоль параллели $x = {{x}_{{i0}}}$ между любыми соседними точками i, j и $i - 1,j$:

Записав решения (1.13), (1.14) между всеми соседними точками в направлениях x и y и мы получим неопределенную систему алгебраических уравнений (СЛАУ). Чтобы полученную СЛАУ сделать определенной, необходимо дополнить ее краевыми условиями и уравнениями связи. Краевые условия устанавливаются следующим образом. Для каждой из точек на контуре x₁₀, x₂₀, …, x_m0 (рис. 1) записываются матричные уравнения вида

где B_x – прямоугольная матрица, выбирающая элементы столбца производных, w_x0j – столбец производных (1.4), v_x0j – столбец задаваемых значений выбранных производных.

Например, для постановки краевого условия защемления, уравнение (1.15) в развернутом виде запишется как:

Для постановки других видов краевых условий следует изменять элементы матрицы B_x. Если же на краю нужно задать не нулевое значение, а, например, приложить силу или изгибающий момент, то следует в столбце v_x0j заменить 0 на необходимое значение. Краевые условия для точек y₀₁, y₀₂, …, y_0n ставятся аналогично.

Дополнив СЛАУ уравнениями связи (1.11), (1.12) для точек пересечения, мы тем получим определенную систему. Для задания поверхностных нагрузок следует использовать одно из уравнений уравнение связи – уравнение (1.11). Заданием в правой части уравнения (1.11) числовых значений величин q₁, q₂, q₃ будут нагружаться соответствующие области вокруг точек пересечения меридианов (рис. 2).

Решив определенную СЛАУ, получаем значения элементов столбцов (1.4), (1.5), а так же значения функций $f_{1}^{{\left( 1 \right)}}(x,y)$ … $f_{3}^{{\left( 3 \right)}}(x,y)$. Решение носит приближенный характер. При сгущении сетки меридианов и параллелей, т.е при увеличении числа точек выполнения уравнений связи и уменьшения областей осреднения функций $f_{1}^{{\left( 1 \right)}}(x,y)$ … $f_{3}^{{\left( 3 \right)}}(x,y)$ решение стремится к аналитическому, то есть с контролируемой погрешностью.

2. Результаты вычислений. Решена краевая задача механики деформирования пологой защемленной по краям оболочки, нагруженной внешним давлением. Расчетная схема представлена на рис. 3. Вычислены перемещения – w – и изгибающие моменты – M – в полюсе пологой оболочки двоякой кривизны. Внешнее давление 10 МПа. Параметры оболочки: размеры в плане – 1 × 1 м, радиусы кривизны – 10 м, толщина – 5 мм. Материал панели – алюминий с характеристиками: модуль Юнга – 0.7 × 10¹¹ Па, коэффициент Пуассона – 0.3. На оболочку накладывалась сетка меридианов и параллелей при разбиении сторон оболочки от 5 × 5 интервалов до 25 × 25 интервалов. Результаты представлены в таблице 1 и в виде графиков. Помимо величин прогибов и моментов (рис. 4, 6) для определения сходимости решения и определения погрешности вычислялась разность их значений полученных на двух соседних итерациях, отнесенная к величине вычисленного значения на данной итерации. В таблице 1 и на графиках (рис. 5, 7) эта величина названа “относительная разность”.

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