Three-Field FEM in Shell Calculations with Options for Interpolation of the Sought Values

Full Text

1. Введение

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

Теория расчета тонких оболочек в настоящее время является достаточно развитой [1–4] при различных видах нагружения. Аналитические решения уравнений теории тонких оболочек оказались возможными лишь в некоторых частных случаях, далеких от инженерной практики. Поэтому актуальными оказались разработки численных методов решения уравнений теории оболочек. Среди численных методов решения задач деформирования инженерных конструкций получил развитие и широкое использование метод конечных элементов (МКЭ) [5–7]. Наиболее широкое использование при расчете оболочек получил МКЭ в формулировке метода перемещений [8–14]. МКЭ в формулировке метода перемещений использовался в расчетах трехслойных оболочек [15–22]. Для обоснования разработанной методики расчета трехслойных оболочек в [21, 22] выполнено сравнение полученных результатов с аналитическими решениями других исследователей.

В расчетах оболочек МКЭ широко использовался и в смешанной формулировке [23–26]. В последние годы разрабатывается виртуальный МКЭ с объемными конечными элементами на пространственных сетках дискретизации [27, 28].

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

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

Получение аппроксимирующих выражений искомых величин с учетом смещения как твердого тела для конечного элемента в форме цилиндрической панели показано в [29, 30], для трехслойных оболочек в [15, 16].

Проблеме учета жестких смещений в МКЭ в векторной формулировке посвящены работы [31, 32].

С целью решения проблем численного анализа НДС оболочек в настоящей работе излагается разработанный трехпольный вариант МКЭ со следующими полями узловых неизвестных:

• перемещения и их первые производные;
• деформации и искривления срединной поверхности;
• усилия и моменты срединной поверхности.

Аппроксимация искомых величин использована в двух вариантах.

В первом варианте традиционные аппроксимирующие выражения использовались для аппроксимации непосредственно компонент векторов перемещений и компонент тензоров деформаций, искривлений, усилий и моментов.

Во втором варианте традиционные аппроксимирующие процедуры применялись к искомым векторам и тензорам, и после координатных преобразований в используемой криволинейной системе координат определялись аппроксимирующие выражения для компонент искомых векторов и тензоров.

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

2. Геометрия оболочки

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

${\overrightarrow{R}}^0=x\overrightarrow{i}+r\left(x, \theta \right)\sin \theta \overrightarrow{j}+r\left(x, \theta \right)\cos \theta \overrightarrow{k}$, (2.1)

где θ – полярный угол, отсчитываемый в плоскости yOz от оси Oz против хода часовой стрелки; r = r (x, θ) – функция, зависящая от конкретного типа поверхности.

Так, например, для оболочки в форме эллипсоида данная функция имеет вид

$r\left(x, \theta \right)=\frac{bc\sqrt{a^2-x^2}}{a\sqrt{b^2{\cos }^2\theta +c^2{\sin }^2\theta }}$, (2.2)

где a, b, c – параметры эллипсоида.

Для эллиптического цилиндра эта формула зависит только от полярного угла

$r\left(\theta \right)=\frac{bc}{\sqrt{b^2{\cos }^2\theta +c^2{\sin }^2\theta }}$ (2.3)

Базисные векторы точки M⁰ срединной поверхности оболочки в исходном состоянии можно получить стандартным образом с использованием формул дифференциальной геометрии

${\overrightarrow{a}}_1^0={\overrightarrow{R}}_x^0, {\overrightarrow{a}}_2^0={\overrightarrow{R}}_{\theta }^0, {\overrightarrow{a}}^0={\overrightarrow{a}}_1^0\times {\overrightarrow{a}}_2^0/\sqrt{a^0}$, (2.4)

где $a^0=\left({\overrightarrow{a}}_1^0\cdot {\overrightarrow{a}}_1^0\right)\left({\overrightarrow{a}}_2^0\cdot {\overrightarrow{a}}_2^0\right)-{\left({\overrightarrow{a}}_1^0\cdot {\overrightarrow{a}}_2^0\right)}^2$.

Векторы базиса точки M^0ζ, расположенной в произвольном слое оболочки на расстоянии ζ от точки M⁰ в исходном состоянии определяются по формулам

${\overrightarrow{g}}_1^0={\overrightarrow{R}}_x^{0\zeta }, {\overrightarrow{g}}_2^0={\overrightarrow{R}}_{\theta }^{0\zeta }$, (2.5)

где ${\overrightarrow{R}}^{0\zeta }={\overrightarrow{R}}^0+\zeta {\overrightarrow{a}}^0$.

В процессе деформирования оболочки под действием приложенной внешней нагрузки точки M⁰ и M^0ζ займут новые положения M и M^ζ, определяемые векторами перемещений $\overrightarrow{v}$ и $\overrightarrow{V}$ соответственно

$\overrightarrow{v}=v^1{\overrightarrow{a}}_1^0+v^2{\overrightarrow{a}}_2^0+v{\overrightarrow{a}}^0, \overrightarrow{V}=\overrightarrow{v}+\zeta \left(\overrightarrow{a}-{\overrightarrow{a}}^0\right)$ (2.6)

Входящий в структуру формулы вектора $\overrightarrow{V}$ орт нормали к деформированной срединной поверхности определяется векторным произведением

$\overrightarrow{a}={\overrightarrow{a}}_1\times {\overrightarrow{a}}_2/\sqrt{a}$, (2.7)

где ${\overrightarrow{a}}_1={\overrightarrow{R}}_x={\left({\overrightarrow{R}}^0+\overrightarrow{v}\right)}_x, {\overrightarrow{a}}_2={\overrightarrow{R}}_{,\theta }={\left({\overrightarrow{R}}^0+\overrightarrow{v}\right)}_{,\theta }, \overrightarrow{a}=\left({\overrightarrow{a}}_1\cdot {\overrightarrow{a}}_1\right)\left({\overrightarrow{a}}_2\cdot {\overrightarrow{a}}_2\right)-{\left({\overrightarrow{a}}_1\cdot {\overrightarrow{a}}_2\right)}^2$.

Векторы базиса в точке M^ζ деформированного состояния оболочки определяются зависимостями

$\begin{array}{l}{\overrightarrow{g}}_1={\overrightarrow{R}}_x^{\zeta }={\left({\overrightarrow{R}}^{0\zeta }+\overrightarrow{V}\right)}_x={\overrightarrow{g}}_1^0+{\overrightarrow{v}}_x+\zeta {\left(\overrightarrow{a}-{\overrightarrow{a}}^0\right)}_x\\ {\overrightarrow{g}}_2={\overrightarrow{R}}_{\theta }^{\zeta }={\left({\overrightarrow{R}}^{0\zeta }+\overrightarrow{V}\right)}_{\theta }={\overrightarrow{g}}_2^0+{\overrightarrow{v}}_{\theta }+\zeta {\left(\overrightarrow{a}-{\overrightarrow{a}}^0\right)}_{\theta }\end{array}$ (2.8)

Метрические тензоры точек M^0ζ и M^ζ определяются компонентами

$g_{\alpha \beta }^0={\overrightarrow{g}}_{\alpha }^0\cdot {\overrightarrow{g}}_{\beta }^0, g_{\alpha \beta }={\overrightarrow{g}}_{\alpha }\cdot {\overrightarrow{g}}_{\beta }$, (2.9)

где нижние индексы α и β последовательно принимают значения 1, 2.

Ковариантные компоненты тензора деформаций в точке M^ζ могут быть получены посредством использования соотношения механики сплошной среды

${\epsilon }_{\alpha \beta }^{\zeta }=\left(g_{\alpha \beta }-g_{\alpha \beta }^0\right)/2$. (2.10)

При использовании гипотезы тонких оболочек Кирхгофа–Лява [1] соотношения (2.10) могут быть представлены суммой

$\underset{3\times 1}{\left\{{\epsilon }_{\alpha \beta }^{\zeta }\right\}}=\underset{3\times 1}{\left\{{\epsilon }_{\alpha \beta }\right\}}+\zeta \underset{3\times 1}{\left\{{\aleph }_{\alpha \beta }\right\}}$, (2.11)

где ${\underset{1\times 3}{\left\{{\epsilon }_{\alpha \beta }^{\zeta }\right\}}}^T=\left\{{\epsilon }_{11}^{\zeta }{\epsilon }_{22}^{\zeta }2{\epsilon }_{12}^{\zeta }\right\}, {\underset{1\times 3}{\left\{{\epsilon }_{\alpha \beta }\right\}}}^T=\left\{{\epsilon }_{11}{\epsilon }_{22}2{\epsilon }_{12}\right\}, {\underset{1\times 3}{\left\{{\aleph }_{\alpha \beta }\right\}}}^T=\left\{{\aleph }_{11}{\aleph }_{22}2{\aleph }_{12}\right\}$ – деформации и искривления в точках M^ζ и M оболочки.

3. Дискретная модель объекта

Для построения дискретной модели оболочки был выбран четырехузловой конечный элемент в виде фрагмента срединной поверхности с узлами i, j, k, l.

Криволинейные координаты x и θ дискретного элемента определяются с использованием билинейных функций через координаты узлов посредством матричных соотношений

$x={\underset{1\times 4}{\left\{\phi \right\}}}^T\underset{4\times 1}{\left\{x_y\right\}}, \theta ={\underset{1\times 4}{\left\{\phi \right\}}}^T\underset{4\times 1}{\left\{{\theta }_y\right\}}$, (3.1)

где ${\underset{1\times 4}{\left\{x_y\right\}}}^T=\left\{x^i x^j x^k x^l\right\}, {\underset{1\times 4}{\left\{{\theta }_y\right\}}}^T=\left\{{\theta }^i {\theta }^j {\theta }^k {\theta }^l\right\}$.

В качестве искомых кинематических и силовых параметров конечного элемента выбираются следующие столбцы узловых неизвестных: ${\underset{1\times 36}{\left\{v_y^L\right\}}}^T=\left\{{\underset{1\times 12}{\left\{v_y^{1L}\right\}}}^T{\underset{1\times 12}{\left\{v_y^{2L}\right\}}}^T{\underset{1\times 12}{\left\{v_y^L\right\}}}^T\right\}$ – строки узловых перемещений, каждая из которых имеет вид ${\underset{1\times 12}{\left\{q_y^L\right\}}}^T=\left\{q^i q^j q^k q^l q_{\xi }^i...q_{\xi }^l q_{\eta }^i\mathrm{...}{q}_{\eta }^l\right\}$, где под q понимается компонента вектора перемещения v¹, v² или v; ${\underset{1\times 24}{\left\{\epsilon {\aleph }_y\right\}}}^T=\left\{{\underset{1\times 12}{\left\{{\epsilon }_y\right\}}}^T{\underset{1\times 12}{\left\{{\aleph }_y\right\}}}^T\right\}$ – строки узловых деформаций и искривлений срединной поверхности ${\underset{1\times 12}{\left\{{\lambda }_y\right\}}}^T=\left\{{\lambda }_{11}^i {\lambda }_{11}^j {\lambda }_{11}^k {\lambda }_{11}^l {\lambda }_{22}^i\mathrm{...}{\lambda }_{22}^l 2{\lambda }_{12}^i\mathrm{...2}{\lambda }_{12}^l\right\}$, где под λ_αβ понимается деформация ε_αβ или искривления ℵ_αβ срединной поверхности; ${\underset{1\times 24}{\left\{N{M}_y\right\}}}^T=\left\{{\underset{1\times 12}{\left\{N_y\right\}}}^T{\underset{1\times 12}{\left\{M_y\right\}}}^T\right\}$ – строки узловых продольных сил и изгибающих моментов ${\underset{1\times 12}{\left\{T_y\right\}}}^T=\left\{T^{11i} T^{11j} T^{11k} T^{11l} T^{22i}\mathrm{...}{T}^{22l} T^{12i}\mathrm{...}{T}^{12l}\right\}$. Здесь под T^αβ понимается контравариантная компонента тензора внутренних усилий или тензора моментов.

В разработанных к настоящему времени вычислительных комплексах (ANSYS, NASTRAN, ABAQUS и других) и конечно-элементных алгоритмах традиционно используется покомпонентная интерполяционная процедура, отличительной особенностью которой является применение интерполяционного выражения отдельно к каждой компоненте вектора перемещения или к каждой компоненте тензорных величин, таких как тензор деформаций, искривлений, сил и моментов. Согласно этой интерполяционной процедуре, можно записать следующие интерполяционные выражения

$q={\underset{1\times 12}{\left\{\psi \right\}}}^T\underset{12\times 1}{\left\{q_y^L\right\}}, {\lambda }_{\alpha \beta }={\underset{1\times 4}{\left\{\phi \right\}}}^T\underset{4\times 1}{\left\{{\lambda }_{\alpha \beta y}\right\}}, T^{\alpha \beta }={\underset{1\times 4}{\left\{\phi \right\}}}^T\underset{4\times 1}{\left\{T_y^{\alpha \beta }\right\}}$, (3.2)

где ${\underset{1\times 12}{\left\{\psi \right\}}}^T$ – матрица-строка, содержащая произведения полиномов Эрмита третьей степени [24, 25].

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

С использованием интерполяционных зависимостей (3.2) формируются аппроксимирующие матричные выражения

$\underset{6\times 1}{\left\{\epsilon \aleph \right\}}=\underset{6\times 24}{\left[H_1\right]}\underset{24\times 1}{\left\{\epsilon {\aleph }_y\right\}}, \underset{6\times 1}{\left\{NM\right\}}=\underset{6\times 24}{\left[H_1\right]}\underset{24\times 1}{\left\{N{M}_y\right\}}, \underset{3\times 1}{\left\{v\right\}}=\underset{3\times 36}{\left[A_1\right]}\underset{36\times 1}{\left\{v_y\right\}}, \underset{6\times 1}{\left\{\epsilon {\aleph }^k\right\}}=\underset{6\times 36}{\left[B_1\right]}\underset{36\times 1}{\left\{v_y\right\}}$, (3.3)

где ${\underset{1\times 6}{\left\{\epsilon \aleph \right\}}}^T=\left\{\underset{1\times 3}{{\left\{{\epsilon }_{\alpha \beta }\right\}}^T}\underset{1\times 3}{{\left\{{\aleph }_{\alpha \beta }\right\}}^T}\right\}$ – строка деформаций и искривлений внутренней точки конечного элемента; ${\underset{1\times 6}{\left\{NM\right\}}}^T=\left\{N^{11} N^{22} N^{12} M^{11} M^{22} M^{12}\right\}$ – строка усилий и моментов внутренней точки конечного элемента; $\underset{1\times 6}{\left\{\epsilon {\aleph }^k\right\}}$ – строка деформаций и искривлений внутренней точки, определяемых функциями перемещений.

Как следует из (3.2), каждая из компонент вектора и тензора интерполируется через узловые значения этой же самой компоненты. Данная аппроксимация является корректной, если векторы и тензоры определены в декартовой системе координат.

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

При реализации трехпольного варианта МКЭ применительно к анализу НДС оболочек тензорно-векторная форма интерполяционной процедуры предусматривает использование следующих интерполяционных выражений

$\tilde{\epsilon }={\underset{1\times 4}{\left\{\phi \right\}}}^T\underset{4\times 1}{\left\{{\tilde{\epsilon }}_y\right\}}, \tilde{\aleph }={\underset{1\times 4}{\left\{\phi \right\}}}^T\underset{4\times 1}{\left\{{\tilde{\aleph }}_y\right\}}, \tilde{N}={\underset{1\times 4}{\left\{\phi \right\}}}^T\underset{4\times 1}{\left\{{\tilde{N}}_y\right\}}, \tilde{M}={\underset{1\times 4}{\left\{\phi \right\}}}^T\underset{4\times 1}{\left\{{\tilde{M}}_y\right\}}, \overrightarrow{v}={\underset{1\times 12}{\left\{\psi \right\}}}^T\underset{12\times 1}{\left\{{\overrightarrow{v}}_y^L\right\}}$, (3.4)

где ${\underset{1\times 4}{\left\{{\tilde{\epsilon }}_y\right\}}}^T=\left\{{\tilde{\epsilon }}^i {\tilde{\epsilon }}^j {\tilde{\epsilon }}^k {\tilde{\epsilon }}^l\right\}; {\underset{1\times 4}{\left\{{\tilde{\aleph }}_y\right\}}}^T=\left\{{\tilde{\aleph }}^i {\tilde{\aleph }}^j {\tilde{\aleph }}^k {\tilde{\aleph }}^l\right\}$ – строки узловых тензоров деформаций и искривлений срединной поверхности; ${\underset{1\times 4}{\left\{{\tilde{N}}_y\right\}}}^T=\left\{{\tilde{N}}^i {\tilde{N}}^j {\tilde{N}}^k {\tilde{N}}^l\right\}$; ${\underset{1\times 4}{\left\{{\tilde{M}}_y\right\}}}^T=\left\{{\tilde{M}}^i {\tilde{M}}^j {\tilde{M}}^k {\tilde{M}}^l\right\}$ – строки узловых тензоров продольных сил и изгибающих моментов; ${\underset{1\times 12}{\left\{{\overrightarrow{v}}_y^L\right\}}}^T=\left\{{\overrightarrow{v}}^i {\overrightarrow{v}}^j {\overrightarrow{v}}^k {\overrightarrow{v}}^l {\overrightarrow{v}}_{\xi }^i\mathrm{...}{\overrightarrow{v}}_{\xi }^l {\overrightarrow{v}}_{\eta }^i\mathrm{...}{\overrightarrow{v}}_{\eta }^l\right\}$ – строка узловых векторов перемещения и их производных первого порядка.

Тензоры деформаций и искривлений, усилий и моментов, входящих в (3.4), можно представить компонентами диадных произведений соответствующих векторов локальных базисов, например

$\begin{array}{l}\tilde{\epsilon }={\epsilon }_{11}{\overrightarrow{a}}^{01}\cdot {\overrightarrow{a}}^{01}+{\epsilon }_{22}{\overrightarrow{a}}^{02}\cdot {\overrightarrow{a}}^{02}+2{\epsilon }_{12}{\overrightarrow{a}}^{01}\cdot {\overrightarrow{a}}^{02}={\underset{1\times 3}{\left\{{\overrightarrow{\overrightarrow{a}}}^{0\alpha \beta }\right\}}}^T\underset{3\times 1}{\left\{{\epsilon }_{\alpha \beta }\right\}}, \tilde{\aleph }={\underset{1\times 3}{\left\{{\overrightarrow{\overrightarrow{a}}}^{0\alpha \beta }\right\}}}^T\underset{3\times 1}{\left\{{\aleph }_{\alpha \beta }\right\}}\\ \tilde{N}=N^{11}{\overrightarrow{a}}_1^0\cdot {\overrightarrow{a}}_1^0+N^{22}{\overrightarrow{a}}_2^0\cdot {\overrightarrow{a}}_2^0+N^{12}{\overrightarrow{a}}_1^0\cdot {\overrightarrow{a}}_2^0={\underset{1\times 3}{\left\{{\overrightarrow{\overrightarrow{a}}}_{\alpha \beta }^0\right\}}}^T\underset{3\times 1}{\left\{N^{\alpha \beta }\right\}}, \tilde{M}={\underset{1\times 3}{\left\{{\overrightarrow{\overrightarrow{a}}}_{\alpha \beta }^0\right\}}}^T\underset{3\times 1}{\left\{M^{\alpha \beta }\right\}}\\ \underset{4\times 1}{\left\{{\tilde{\epsilon }}_y\right\}}=\underset{4\times 12}{\left[\overrightarrow{\overrightarrow{D}}\right]}\underset{12\times 1}{\left\{{\epsilon }_{\alpha \beta y}\right\}}, \underset{4\times 1}{\left\{{\tilde{\aleph }}_y\right\}}=\underset{4\times 12}{\left[\overrightarrow{\overrightarrow{D}}\right]}\underset{12\times 1}{\left\{{\aleph }_{\alpha \beta y}\right\}}, \underset{4\times 1}{\left\{{\tilde{N}}_y\right\}}=\underset{4\times 12}{\left[\overrightarrow{\overrightarrow{L}}\right]}\underset{12\times 1}{\left\{N_y^{\alpha \beta }\right\}}, \underset{4\times 1}{\left\{{\tilde{M}}_y\right\}}=\underset{4\times 12}{\left[\overrightarrow{\overrightarrow{L}}\right]}\underset{12\times 1}{\left\{M_y^{\alpha \beta }\right\}},\end{array}$ (3.5)

где

$\underset{4\times 12}{\left[\overrightarrow{\overrightarrow{L}}\right]}=\left[\begin{array}{cccccccccccc}{\left({\overrightarrow{a}}_1^0{\overrightarrow{a}}_1^0\right)}^i& {\left({\overrightarrow{a}}_2^0{\overrightarrow{a}}_2^0\right)}^i& {\left({\overrightarrow{a}}_1^0{\overrightarrow{a}}_2^0\right)}^i& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& {\left({\overrightarrow{a}}_1^0{\overrightarrow{a}}_1^0\right)}^j& {\left({\overrightarrow{a}}_2^0{\overrightarrow{a}}_2^0\right)}^j& {\left({\overrightarrow{a}}_1^0{\overrightarrow{a}}_2^0\right)}^j& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& {\left({\overrightarrow{a}}_1^0{\overrightarrow{a}}_1^0\right)}^k& {\left({\overrightarrow{a}}_2^0{\overrightarrow{a}}_2^0\right)}^k& {\left({\overrightarrow{a}}_1^0{\overrightarrow{a}}_2^0\right)}^k& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& {\left({\overrightarrow{a}}_1^0{\overrightarrow{a}}_1^0\right)}^l& {\left({\overrightarrow{a}}_2^0{\overrightarrow{a}}_2^0\right)}^l& {\left({\overrightarrow{a}}_1^0{\overrightarrow{a}}_2^0\right)}^l\end{array}\right]$

матрица $\underset{4\times 12}{\left[\overrightarrow{\overrightarrow{D}}\right]}$ имеет вид $\underset{4\times 12}{\left[\overrightarrow{\overrightarrow{L}}\right]}$, но ее компоненты, отличные от нуля, являются диадными произведениями контравариантных векторов базисов узловых точек;

$$\begin{gathered} {\underset{1\times 12}{\left\{{\epsilon }_{\alpha \beta y}\right\}}}^T=\left\{\underset{1\times 3}{\left\{{\epsilon }_{11}^i {\epsilon }_{22}^i 2{\epsilon }_{12}^i\right\}}\underset{1\times 3}{\left\{{\epsilon }_{11}^j {\epsilon }_{22}^j 2{\epsilon }_{12}^j\right\}}\underset{1\times 3}{\left\{{\epsilon }_{11}^k {\epsilon }_{22}^k 2{\epsilon }_{12}^k\right\}}\underset{1\times 3}{\left\{{\epsilon }_{11}^l {\epsilon }_{22}^l 2{\epsilon }_{12}^l\right\}}\right\} \\ {\underset{1\times 12}{\left\{{\aleph }_{\alpha \beta y}\right\}}}^T=\left\{\underset{1\times 3}{\left\{{\aleph }_{11}^i {\aleph }_{22}^i 2{\aleph }_{12}^i\right\}}\underset{1\times 3}{\left\{{\aleph }_{11}^j {\aleph }_{22}^j 2{\aleph }_{12}^j\right\}}\underset{1\times 3}{\left\{{\aleph }_{11}^k {\aleph }_{22}^k 2{\aleph }_{12}^k\right\}}\underset{1\times 3}{\left\{{\aleph }_{11}^l {\aleph }_{22}^l 2{\aleph }_{12}^l\right\}}\right\} \\ {\underset{1\times 12}{\left\{N_y^{\alpha \beta }\right\}}}^T=\left\{\underset{1\times 3}{\left\{N^{11i} N^{22i} N^{12i}\right\}}\underset{1\times 3}{\left\{N^{11j} N^{22j} N^{12j}\right\}}\underset{1\times 3}{\left\{N^{11k} N^{22k} N^{12k}\right\}}\underset{1\times 3}{\left\{N^{11l} N^{22l} N^{12l}\right\}}\right\} \\ {\underset{1\times 12}{\left\{M_y^{\alpha \beta }\right\}}}^T=\left\{\underset{1\times 3}{\left\{M^{11i} M^{22i} M^{12i}\right\}}\underset{1\times 3}{\left\{M^{11j} M^{22j} M^{12j}\right\}}\underset{1\times 3}{\left\{M^{11k} M^{22k} M^{12k}\right\}}\underset{1\times 3}{\left\{M^{11l} M^{22l} M^{12l}\right\}}\right\} \end{gathered}$$

Векторы базисов внутренних и узловых точек конечного элемента выражаются через орты декартовой системы координат посредством следующих матричных соотношений

$\underset{3\times 1}{\left\{{\overrightarrow{a}}_{\alpha }^0\right\}}=\underset{3\times 3}{\left[m^0\right]}\underset{3\times 1}{\left\{\overrightarrow{i}\right\}}, \underset{3\times 1}{\left\{{\overrightarrow{a}}_{\alpha }^{0\rho }\right\}}=\underset{3\times 3}{\left[m^{0\rho }\right]}\underset{3\times 1}{\left\{\overrightarrow{i}\right\}}$, (3.6)

где ${\underset{1\times 3}{\left\{{\overrightarrow{a}}_{\alpha }^0\right\}}}^T=\left\{{\overrightarrow{a}}_1^0 {\overrightarrow{a}}_1^0 {\overrightarrow{a}}^0\right\}$, ${\underset{1\times 3}{\left\{\overrightarrow{i}\right\}}}^T=\left\{\overrightarrow{i} \overrightarrow{j} \overrightarrow{k}\right\}$, ${\underset{1\times 3}{\left\{{\overrightarrow{a}}_{\alpha }^{0\rho }\right\}}}^T=\left\{{\overrightarrow{a}}_1^{0\rho } {\overrightarrow{a}}_2^{0\rho } {\overrightarrow{a}}^{0\rho }\right\}$. Верхний индекс ρ обозначает конкретный из узлов КЭ i, j, k, l.

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

$\underset{3\times 1}{\left\{{\overrightarrow{a}}_{\alpha }^{0\rho }\right\}}=\underset{3\times 3}{\left[m^{0\rho }\right]}{\underset{3\times 3}{\left[m^0\right]}}^{-1}\underset{3\times 1}{\left\{{\overrightarrow{a}}_{\alpha }^0\right\}}=\underset{3\times 3}{\left[{\omega }^{\rho }\right]}\underset{3\times 1}{\left\{{\overrightarrow{a}}_{\alpha }^0\right\}}$ (3.7)

С учетом (3.7) диадные произведения узлов КЭ также могут быть выражены через диадные произведения точки, принадлежащей внутренней области конечного элемента

$\underset{3\times 1}{\left\{{\overrightarrow{\overrightarrow{a}}}_{\alpha \beta }^{0\rho }\right\}}=\underset{3\times 3}{\left[b^{0\rho }\right]}\underset{3\times 1}{\left\{{\overrightarrow{\overrightarrow{a}}}_{\alpha \beta }^0\right\}}, \underset{3\times 1}{\left\{{\overrightarrow{\overrightarrow{a}}}^{0\alpha \beta \rho }\right\}}=\underset{3\times 3}{\left[d^{0\rho }\right]}\underset{3\times 1}{\left\{{\overrightarrow{\overrightarrow{a}}}^{0\alpha \beta }\right\}}$, (3.8)

где ${\underset{1\times 3}{\left\{{\overrightarrow{\overrightarrow{a}}}_{\alpha \beta }^{0\rho }\right\}}}^T=\left\{\left({\overrightarrow{a}}_1^{0\rho }{\overrightarrow{a}}_1^{0\rho }\right)\left({\overrightarrow{a}}_2^{0\rho }{\overrightarrow{a}}_2^{0\rho }\right)\left({\overrightarrow{a}}_1^{0\rho }{\overrightarrow{a}}_2^{0\rho }\right)\right\}, {\underset{1\times 3}{\left\{{\overrightarrow{\overrightarrow{a}}}^{0\alpha \beta \rho }\right\}}}^T=\left\{\left({\overrightarrow{a}}^{01}{\overrightarrow{a}}^{01}\right)\left({\overrightarrow{a}}^{02}{\overrightarrow{a}}^{02}\right)\left({\overrightarrow{a}}^{01}{\overrightarrow{a}}^{02}\right)\right\}$.

Принимая во внимание (3.8), матрица $\left[\overrightarrow{\overrightarrow{L}}\right]$ примет вид

$\underset{4\times 12}{\left[\overrightarrow{\overrightarrow{L}}\right]}=\left[\begin{array}{ccccccccccc}{\left\{{\overrightarrow{\overrightarrow{a}}}_{\alpha \beta }^0\right\}}^T& {\left[b^{0i}\right]}^T& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& {\left\{{\overrightarrow{\overrightarrow{a}}}_{\alpha \beta }^0\right\}}^T& {\left[b^{0j}\right]}^T& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& {\left\{{\overrightarrow{\overrightarrow{a}}}_{\alpha \beta }^0\right\}}^T& {\left[b^{0k}\right]}^T& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& {\left\{{\overrightarrow{\overrightarrow{a}}}_{\alpha \beta }^0\right\}}^T& {\left[b^{0l}\right]}^T\end{array}\right]$ (3.9)

При учете (3.5)–(3.9) соотношения (3.4) могут быть представлены выражениями

$\begin{array}{l}\underset{1\times 3}{\overset{T}{\left\{{\overrightarrow{\overrightarrow{a}}}^{0\alpha \beta }\right\}}}\underset{3\times 1}{\left\{{\epsilon }_{\alpha \beta }\right\}}=\underset{1\times 3}{\overset{T}{\left\{{\overrightarrow{\overrightarrow{a}}}^{0\alpha \beta }\right\}}}\left[\left.\underset{3\times 3}{{\varphi }_1{\left[d^{0i}\right]}^T}\right|\left.\underset{3\times 3}{{\varphi }_2{\left[d^{0j}\right]}^T}\right|\left.\underset{3\times 3}{{\varphi }_3{\left[d^{0k}\right]}^T}\right|\underset{3\times 3}{{\varphi }_4{\left[d^{0l}\right]}^T}\right]\underset{12\times 1}{\left\{{\epsilon }_{\alpha \beta y}\right\}}\\ \underset{1\times 3}{\overset{T}{\left\{{\overrightarrow{\overrightarrow{a}}}^{0\alpha \beta }\right\}}}\underset{3\times 1}{\left\{{\aleph }_{\alpha \beta }\right\}}=\underset{1\times 3}{\overset{T}{\left\{{\overrightarrow{\overrightarrow{a}}}^{0\alpha \beta }\right\}}}\left[\left.\underset{3\times 3}{{\varphi }_1{\left[d^{0i}\right]}^T}\right|\left.\underset{3\times 3}{{\varphi }_2{\left[d^{0j}\right]}^T}\right|\left.\underset{3\times 3}{{\varphi }_3{\left[d^{0k}\right]}^T}\right|\underset{3\times 3}{{\varphi }_4{\left[d^{0l}\right]}^T}\right]\underset{12\times 1}{\left\{{\aleph }_{\alpha \beta y}\right\}}\\ \underset{1\times 3}{\overset{T}{\left\{{\overrightarrow{\overrightarrow{a}}}_{\alpha \beta }^0\right\}}}\underset{3\times 1}{\left\{N^{\alpha \beta }\right\}}=\underset{1\times 3}{\overset{T}{\left\{{\overrightarrow{\overrightarrow{a}}}_{\alpha \beta }^0\right\}}}\left[\left.\underset{3\times 3}{{\varphi }_1{\left[b^{0i}\right]}^T}\right|\left.\underset{3\times 3}{{\varphi }_2{\left[b^{0j}\right]}^T}\right|\left.\underset{3\times 3}{{\varphi }_3{\left[b^{0k}\right]}^T}\right|\underset{3\times 3}{{\varphi }_4{\left[b^{0l}\right]}^T}\right]\underset{12\times 1}{\left\{N_y^{\alpha \beta }\right\}}\\ \underset{1\times 3}{\overset{T}{\left\{{\overrightarrow{\overrightarrow{a}}}_{\alpha \beta }^0\right\}}}\underset{3\times 1}{\left\{M^{\alpha \beta }\right\}}=\underset{1\times 3}{\overset{T}{\left\{{\overrightarrow{\overrightarrow{a}}}_{\alpha \beta }^0\right\}}}\left[\left.\underset{3\times 3}{{\varphi }_1{\left[b^{0i}\right]}^T}\right|\left.\underset{3\times 3}{{\varphi }_2{\left[b^{0j}\right]}^T}\right|\left.\underset{3\times 3}{{\varphi }_3{\left[b^{0k}\right]}^T}\right|\underset{3\times 3}{{\varphi }_4{\left[b^{0l}\right]}^T}\right]\underset{12\times 1}{\left\{M_y^{\alpha \beta }\right\}}\end{array}$ (3.10)

Из (3.10) можно получить необходимые интерполяционные зависимости для компонент тензоров деформаций, искривлений, продольных сил и изгибающих моментов, например

$\begin{array}{l}\left\{\begin{array}{l}{\epsilon }_{11}\\ {\epsilon }_{22}\\ 2{\epsilon }_{12}\end{array}\right\}=\left\{\left.{\varphi }_1\underset{3\times 3}{\overset{T}{\left[d^{0i}\right]}}\right|\left.{\varphi }_2\underset{3\times 3}{\overset{T}{\left[d^{0j}\right]}}\right|\left.{\varphi }_3\underset{3\times 3}{\overset{T}{\left[d^{0k}\right]}}\right|{\varphi }_4\underset{3\times 3}{\overset{T}{\left[d^{0l}\right]}}\right\}\underset{12\times 1}{\left\{{\epsilon }_{\alpha \beta y}\right\}}\\ \left\{\begin{array}{l}{\aleph }_{11}\\ {\aleph }_{22}\\ 2{\aleph }_{12}\end{array}\right\}=\left\{\left.{\varphi }_1\underset{3\times 3}{\overset{T}{\left[d^{0i}\right]}}\right|\left.{\varphi }_2\underset{3\times 3}{\overset{T}{\left[d^{0j}\right]}}\right|\left.{\varphi }_3\underset{3\times 3}{\overset{T}{\left[d^{0k}\right]}}\right|{\varphi }_4\underset{3\times 3}{\overset{T}{\left[d^{0l}\right]}}\right\}\underset{12\times 1}{\left\{{\aleph }_{\alpha \beta y}\right\}}\end{array}$

$\begin{array}{l}\left\{\begin{array}{l}{N}^{11}\\ N^{22}\\ N^{12}\end{array}\right\}=\left\{\left.{\varphi }_1\underset{3\times 3}{\overset{T}{\left[b^{0i}\right]}}\right|\left.{\varphi }_2\underset{3\times 3}{\overset{T}{\left[b^{0j}\right]}}\right|\left.{\varphi }_3\underset{3\times 3}{\overset{T}{\left[b^{0k}\right]}}\right|{\varphi }_4\underset{3\times 3}{\overset{T}{\left[b^{0l}\right]}}\right\}\underset{12\times 1}{\left\{N_y^{\alpha \beta }\right\}}\\ \left\{\begin{array}{l}{M}^{11}\\ M^{22}\\ M^{12}\end{array}\right\}=\left\{\left.{\varphi }_1\underset{3\times 3}{\overset{T}{\left[b^{0i}\right]}}\right|\left.{\varphi }_2\underset{3\times 3}{\overset{T}{\left[b^{0j}\right]}}\right|\left.{\varphi }_3\underset{3\times 3}{\overset{T}{\left[b^{0k}\right]}}\right|{\varphi }_4\underset{3\times 3}{\overset{T}{\left[b^{0l}\right]}}\right\}\underset{12\times 1}{\left\{M_y^{\alpha \beta }\right\}}\end{array}$ (3.11)

Получение интерполяционных соотношений для компонент вектора перемещения в соответствие с векторной формой интерполяционной процедуры изложено [24, 25].

На основе интерполяционных зависимостей (3.11) формируются аппроксимирующие выражения

$\underset{6\times 1}{\left\{\epsilon \aleph \right\}}=\underset{6\times 24}{\left[H_2\right]}\underset{24\times 1}{\left\{\epsilon {\aleph }_y\right\}}, \underset{6\times 1}{\left\{NM\right\}}=\underset{6\times 24}{\left[H_2\right]}\underset{24\times 1}{\left\{N{M}_y\right\}}, \underset{3\times 1}{\left\{v\right\}}=\underset{3\times 36}{\left[A_2\right]}\underset{36\times 1}{\left\{v_y\right\}}, \underset{6\times 1}{\left\{\xi {\aleph }^k\right\}}=\underset{6\times 36}{\left[B_2\right]}\underset{36\times 1}{\left\{v_y\right\}}$. (3.12)

Анализируя (3.11), можно констатировать, что отдельная компонента любого тензора (деформаций, искривлений, усилий и моментов) в точке, принадлежащей внутренней области используемого элемента дискретизации, зависит от узловых значений всех его компонент, то есть от полного набора узловых варьируемых параметров соответствующего тензора.

При общепринятой форме интерполяционной процедуры (3.2) отдельная компонента тензора, например, компонента тензора искривлений ℵ₁₁ является функцией узловых значений только этой же самой компоненты, т.е. ${\aleph }_{11}=f\left({\aleph }_{11}^i {\aleph }_{11}^j {\aleph }_{11}^k {\aleph }_{11}^l\right)$ и не зависит от узловых значений остальных компонент данного тензора ${\aleph }_{22}^{\rho }, {\aleph }_{12}^{\rho }$.

Также особенностью новых интерполяционных зависимостей (3.11) является присутствие в них параметров применяемой в конкретном расчете оболочки криволинейной системы координат, что не наблюдается в традиционных интерполяционных соотношениях (3.2).

4. Матрица жесткости конечного элемента

Для получения матрицы жесткости конечного элемента используется функционал [24, 25], основанный на равенстве работ заданных сил на перемещениях и внутренних усилий на деформациях и искривлениях, дополненный условием равенства нулю работы на деформациях и искривлениях невязки, определяемой как разность между внутренними усилиями, вводимыми по определению, и внутренними усилиями, выражаемыми через деформации и искривления

${\prod }_L=\underset{F}{\int }\underset{1\times 6}{\overset{T}{\left\{NM\right\}}}\underset{6\times 1}{\left\{\epsilon \aleph \right\}}dF-\underset{F}{\int }\underset{1\times 6}{\overset{T}{\left\{\epsilon \aleph \right\}}}\left[\underset{6\times 1}{\left\{NM\right\}}-\underset{6\times 6}{\left[C\right]}\underset{6\times 1}{\left\{\epsilon {\aleph }^k\right\}}\right]dF-\underset{F}{\int }\underset{1\times 3}{\overset{T}{\left\{U\right\}}}\underset{3\times 1}{\left\{P\right\}}dF$, (4.1)

где $\underset{1\times 6}{\overset{T}{\left\{NM\right\}}}=\left\{N^{11} N^{22} N^{12} M^{11} M^{22} M^{12}\right\}$, $\underset{1\times 6}{\overset{T}{\left\{\epsilon \aleph \right\}}}=\left\{{\epsilon }_{11} {\epsilon }_{22} 2{\epsilon }_{12} {\aleph }_{11} {\aleph }_{22} 2{\aleph }_{12}\right\}$ – матрицы-строки продольных сил и изгибающих моментов, а также деформационных параметров в точке срединной поверхности оболочки; {εℵ^k} – столбец, элементы которого определяются соотношениями Коши; ${\underset{1\times 3}{\left\{U\right\}}}^T=\left\{v^1 v^2 v\right\}$; ${\underset{1\times 3}{\left\{P\right\}}}^T=\left\{p_1 p_2 p_3\right\}$ – матрицы-строки компонент вектора перемещения и вектора внешней поверхностной нагрузки.

Входящая в (4.1) матрица [C] определяет связь между столбцами {NM} и {εℵ}

$\underset{6\times 1}{\left\{NM\right\}}=\underset{6\times 6}{\left[C\right]}\underset{6\times 1}{\left\{\epsilon \aleph \right\}}$ (4.2)

При учете аппроксимирующих выражений для кинематических {U}, {εℵ} и силовых {NM} искомых неизвестных в зависимости от используемого варианта аппроксимации (3.2) или (3.11) функционал (4.1) с учетом (4.2) может быть преобразован к виду

$$\begin{gathered} \prod _L=\underset{1\times 24}{\overset{T}{\left\{N{M}_y\right\}}}\underset{F}{\int }\underset{24\times 6}{\overset{T}{\left[H_{\gamma }\right]}}\underset{6\times 6}{\overset{-1}{\left[C\right]}}\underset{6\times 24}{\left[H_{\gamma }\right]}dF\underset{24\times 1}{\left\{N{M}_y\right\}}-\underset{1\times 24}{\overset{T}{\left\{\epsilon {\aleph }_y\right\}}}\underset{F}{\int }\underset{24\times 6}{\overset{T}{\left[H_{\gamma }\right]}}\underset{6\times 24}{\left[H_{\gamma }\right]}dF\underset{24\times 1}{\left\{N{M}_y\right\}}+ \\ +\underset{1\times 24}{\overset{T}{\left\{\epsilon {\aleph }_y\right\}}}\underset{F}{\int }\underset{24\times 6}{\overset{T}{\left[H_{\gamma }\right]}}\underset{6\times 6}{\left[C\right]}\underset{6\times 36}{\left[B_{\gamma }\right]}dF\underset{36\times 36}{\left[P_R\right]}\underset{36\times 1}{\left\{v_y^G\right\}}-\underset{1\times 36}{\overset{T}{\left\{v_y^G\right\}}}\underset{36\times 36}{\overset{T}{\left[P_R\right]}}\underset{F}{\int }\underset{36\times 3}{\left[A_{\gamma }\right]}\underset{3\times 1}{\left\{P\right\}}dF, \end{gathered}$$ (4.3)

где матрицы $\underset{6\times 24}{\left[H_{\gamma }\right]}$, $\underset{6\times 36}{\left[B_{\gamma }\right]}$ и $\underset{36\times 3}{\left[A_{\gamma }\right]}$ компонуются на основе интерполяционных выражений (3.3) или (3.12), нижний индекс γ принимает значения 1, 2, причем γ = 1 соответствует использованию общепринятой покомпонентной интерполяционной процедуре (3.2), а γ = 2 соответствует применению разработанной тензорно-векторной форме интерполяции искомых величин (3.11) в трехпольном варианте МКЭ.

Входящая в (4.3) матрица $\underset{36\times 36}{\left[P_R\right]}$ определяет связь между столбцами кинематических узловых неизвестных в локальной и глобальной системах координат.

Осуществляя минимизацию (4.3) по кинематическим ${\left\{v_y^G\right\}}^T$, {εℵ_y}^T и силовым {NM_y}^T искомым неизвестным, можно записать систему матричных уравнений

$\begin{array}{l}\frac{\partial {\displaystyle \prod _L}}{\partial {\left\{v_y^G\right\}}^T}=\underset{36\times 36}{\overset{T}{\left[P_R\right]}}\underset{F}{\int }\underset{36\times 6}{\overset{T}{\left[B\right]}}\underset{6\times 6}{\overset{T}{\left[C\right]}}\underset{6\times 24}{\left[H_{\gamma }\right]}dF\underset{24\times 1}{\left\{\epsilon {\aleph }_y\right\}}-\underset{36\times 36}{\overset{T}{\left[P_R\right]}}\underset{F}{\int }\underset{36\times 3}{\left[A_{\gamma }\right]}\underset{3\times 1}{\left\{P\right\}}dF=0\\ \frac{\partial {\displaystyle \prod _L}}{\partial {\left\{\epsilon {\aleph }_y\right\}}^T}=\underset{F}{\int }\underset{24\times 6}{\overset{T}{\left[H_{\gamma }\right]}}\underset{6\times 6}{\left[C\right]}\underset{6\times 36}{\left[B\right]}dF\underset{36\times 36}{\left[P_R\right]}\underset{36\times 1}{\left\{v_y^G\right\}}-\underset{F}{\int }\underset{24\times 6}{\overset{T}{\left[H_{\gamma }\right]}}\underset{6\times 24}{\left[H_{\gamma }\right]}dF\underset{24\times 1}{\left\{N{M}_y\right\}}=0\\ \frac{\partial {\displaystyle \prod _L}}{\partial {\left\{N{M}_y\right\}}^T}=\underset{F}{\int }\underset{24\times 6}{\overset{T}{\left[H_{\gamma }\right]}}\underset{6\times 6}{\overset{-1}{\left[C\right]}}\underset{6\times 24}{\left[H_{\gamma }\right]}dF\underset{24\times 1}{\left\{N{M}_y\right\}}-\underset{F}{\int }\underset{24\times 6}{\overset{T}{\left[H_{\gamma }\right]}}\underset{6\times 24}{\left[H_{\gamma }\right]}dF\underset{24\times 1}{\left\{\epsilon {\aleph }_y\right\}}=0\end{array}$ (4.4)

или в более компактном виде

$\begin{array}{l}\underset{36\times 24}{\overset{T}{\left[S\right]}}\underset{24\times 1}{\left\{\epsilon {\aleph }_y\right\}}-\underset{36\times 1}{\left\{f\right\}}=0\\ \underset{24\times 36}{\left[S\right]}\underset{36\times 1}{\left\{v_y^G\right\}}-\underset{24\times 24}{\left[G\right]}\underset{24\times 1}{\left\{N{M}_y\right\}}=0\\ \underset{24\times 24}{\left[W\right]}\underset{24\times 1}{\left\{N{M}_y\right\}}-\underset{24\times 24}{\overset{T}{\left[G\right]}}\underset{24\times 1}{\left\{\epsilon {\aleph }_y\right\}}=\mathrm{0,}\end{array}$ (4.5)

где $\underset{24\times 36}{\left[S\right]}=\underset{F}{\int }\underset{24\times 6}{\overset{T}{\left[H_{\gamma }\right]}}\underset{6\times 6}{\left[C\right]}\underset{6\times 36}{\left[B\right]}dF\underset{36\times 36}{\left[P_R\right]}, \underset{24\times 24}{\left[G\right]}=\underset{F}{\int }\underset{24\times 6}{\overset{T}{\left[H_{\gamma }\right]}}\underset{6\times 24}{\left[H_{\gamma }\right]}dF, \underset{24\times 24}{\left[W\right]}=\underset{F}{\int }\underset{24\times 6}{\overset{T}{\left[H_{\gamma }\right]}}\underset{6\times 6}{\overset{-1}{\left[C\right]}}\underset{6\times 24}{\left[H_{\gamma }\right]}dF$.

Выражая из третьего и второго уравнений (4.5) столбцы {εℵ_y}^T и {NM_y}^T, можно получить следующие матричные зависимости

$$\begin{gathered} \underset{24\times 1}{\left\{\epsilon {\aleph }_y\right\}}=\underset{24\times 24}{\overset{-1}{\left[{\left[G\right]}^T\right]}}\underset{24\times 24}{\left[W\right]}\underset{24\times 1}{\left\{N{M}_y\right\}}, \underset{24\times 1}{\left\{N{M}_y\right\}}=\underset{24\times 24}{\overset{-1}{\left[G\right]}}\underset{24\times 36}{\left[S\right]}\underset{36\times 1}{\left\{v_y^G\right\}} \\ \underset{24\times 1}{\left\{\epsilon {\aleph }_y\right\}}=\underset{24\times 24}{\overset{-1}{\left[{\left[G\right]}^T\right]}}\underset{24\times 24}{\left[W\right]}\underset{24\times 24}{\overset{-1}{\left[G\right]}}\underset{24\times 36}{\left[S\right]}\underset{36\times 1}{\left\{v_y^G\right\}} \end{gathered}$$ (4.6)

Осуществляя подстановку (4.6) в первое уравнение (4.5), можно записать следующее матричное выражение

$\underset{36\times 36}{\left[K\right]}\underset{36\times 1}{\left\{v_y^G\right\}}=\underset{36\times 1}{\left\{f\right\}}$, (4.7)

где $\underset{36\times 36}{\left[K\right]}={\underset{36\times 24}{\left[S\right]}}^T{\underset{24\times 24}{\left[{\left[G\right]}^{-1}\right]}}^T\underset{24\times 24}{\left[W\right]}{\underset{24\times 24}{\left[G\right]}}^{-1}\underset{24\times 36}{\left[S\right]}$ – искомая матрица жесткости размером 36 × 36 конечного элемента, используемого для построения дискретной модели оболочки.

Глобальная матрица жесткости всей исследуемой оболочки, представляющая собой матрицу системы линейных алгебраических уравнений (СЛАУ), компонуется из отдельных матриц жесткостей [K] посредством использования матрицы индексов [5].

Полученные в результате решения методом Гаусса глобальной СЛАУ столбцы $\left\{v_y^G\right\}$ посредством (4.6), позволяют вычислять требуемые для анализа НДС оболочки столбцы продольных сил и изгибающих моментов {εℵ_y}, а также столбцы деформаций и искривлений {NM_y} в любой интересующей исследователя узловой точке рассчитываемой конструкции.

5. Вычислительные эксперименты

Эксперимент 1. В качестве тестовой задачи был выполнен расчет двухшарнирной арки эллиптического очертания (рис. 1), загруженной линейной нагрузкой интенсивности q = 0.1 Н/см.

Рис. 1. Расчетная схема эллиптического кольца при шарнирном опирании

По условиям симметрии рассматривалась половина арки. Использовались следующие исходные данные: L = 1 см; a = 50 см; b = 10 см; h = 0.1 см; E = 2 ×10⁵ МПа;
ν = 0.3. Расчеты выполнялись в двух вариантах: в первом варианте искомые величины аппроксимировались как составляющие скалярных полей; во втором варианте интерполяционные функции применялись к векторным и тензорным полям и только после координатных преобразований определялись аппроксимирующие выражения искомых величин. Результаты повариантных расчетов приведены в табл. 1, где в зависимости от степени сгущения сетки дискретизации даются значения нормальных напряжений на внутренней σⁱⁿ, наружной σ^out и срединной σ^midl поверхностях оболочки, а также моментов и продольных сил в точках арки. В правой колонке табл. 1 даются значения $M_{ф}^{22}$ и $N_{ф}^{22}$, определенные соотношениями

$$\begin{gathered} N_{ф}^{22}=-q\cdot L=-0.1\cdot 1=-0.100 Н \\ {\sigma }_{\theta \theta }^{midl}=-N_{ф}^{22}/\left(L\cdot h\right)=-0.100 Н/\left(1 см\cdot 0.1 см\right)=-1.000 Н/с{м}^2 \\ {M}_{ф}^{22}=0.000 Н\cdot см \end{gathered}$$

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

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

При замене шарнирной опоры на пружинную (рис. 2) оболочка получает возможность смещаться в вертикальном направлении при неизменном НДС.

Рис. 2. Расчетная схема эллиптического кольца при пружинном опирании

В табл. 2 приведены значения параметров НДС оболочки при различных смещениях арки по вертикали. При выполнении повариантных расчетов была использована сетка дискретизации 211 × 2.

Таблица 2. Значения нормальных напряжений, изгибающих моментов и продольных сил в точке приложения нагрузки и в точке пружинного опирания

Анализ данных табл. 2 показывает, что в первом варианте расчета параметры НДС оболочки в точке пружинного опирания претерпевают весьма значительные изменения, которые возрастают по мере увеличения смещения оболочки как абсолютно твердого тела. Так, при смещении оболочки на величину 0.10 м напряжения ${\sigma }_{\theta \theta }^{in}$ уже меняют свой знак, а σ_xx возрастают на несколько порядков. На основе этого можно отметить, что использование общепринятой в МКЭ интерполяции компонент вектора перемещения, компонент тензоров усилий и моментов, деформаций и искривлений не позволяет учитывать смещения оболочки как абсолютно твердого тела. И напротив, использование разработанной тензорно-векторной формы интерполяционной процедуры в трехпольной формулировке МКЭ позволяет в полной мере учитывать смещения оболочки как абсолютно твердого тела. Подтверждением этого факта является анализ значений параметров НДС оболочки, полученных во втором варианте расчета. Анализ данных, представленных в правой половине таблицы 2 показывает, что численные значения нормальных напряжений, продольных сил и изгибающих моментов остаются неизменными даже при существенной величине жестких смещений, что доказывает несомненные преимущества разработанного трехпольного конечного элемента с тензорно-векторной формой интерполяционной процедуры.

Эксперимент 2. Определено НДС эллиптической цилиндрической оболочки, имеющей первоначально шарнирное опирание по торцам (рис. 3).

Рис. 3. Расчетная схема эллиптического цилиндра при шарнирном опирании

Вдоль верхней образующей была приложена линейно распределенная нагрузка интенсивности q = 10 Н/см. Геометрические размеры оболочки были приняты равными: L = 2.0 м; толщина стенки h = 0.015 м; параметры эллипса поперечного сечения: b = 1.0 м; c = 0.3 м. Физические характеристики: E = 2 × 10⁵ МПа; ν = 0.3.

Как и в эксперименте 1 определение НДС оболочки выполнялось при двух вариантах аппроксимации искомых величин. Результаты повариантных расчетов при шарнирном способе опирания оболочки приведены в табл. 3, в которой даются численные значения физических напряжений σ_θθ и σ_xx на внутренней и наружной поверхностях оболочки в точках A, B, C и D (рис. 3) при различных размерах сетки дискретизации ¼ части эллиптического цилиндра (была рассчитана только четвертая часть оболочки вследствие симметрии расчетной схемы).

Таблица 3. Значения нормальных напряжений в точках при шарнирном опирании

Анализ данных, приведенных в табл. 3, показывает, что при использовании обоих вариантов аппроксимации искомых величин наблюдается сходимость вычислительного процесса, но во втором варианте аппроксимации искомых величин сходимость выше. Значения параметров НДС в обоих вариантах весьма близки между собой.

При замене шарнирных опор оболочки на пружинные при действии той же нагрузки она получает возможность смещаться по вертикали ка твердое тело (рис. 4).

Рис. 4. Расчетная схема эллиптического цилиндра при пружинном опирании

Результаты определения НДС оболочки при фиксированной сетке узлов  и различных смещениях оболочки как твердого тела приведены в табл. 4.

Таблица 4. Значения нормальных напряжений в точках при пружинном опирании