Расчет оболочек вращения при использовании смешанного МКЭ с векторной аппроксимационной процедурой

Толық мәтін

Анализ напряженно-деформированного состояния (НДС) объектов техносферы в виде оболочек вращения и их фрагментов в настоящее время осуществляют с помощью численных методов [1, 2], доминирующее положение среди которых занимает метод конечных элементов (МКЭ) [3–5]. Подавляющее большинство современных зарубежных и отечественных вычислительных комплексов типа ANSYS, NASTRAN, ABAQUS, “Лира” реализуют конечно-элементную процедуру в форме метода перемещений. В то же время использование МКЭ в смешанной формулировке [6] позволяет создавать вычислительные алгоритмы, обладающие рядом преимуществ по сравнению с алгоритмами, построенными на базе метода перемещений [7, 8]. В настоящей статье реализован конечно-элементный алгоритм анализа НДС оболочек вращения при реализации смешанного варианта МКЭ с использованием векторного способа интерполяции кинематических искомых неизвестных. На тестовых примерах оболочки вращения, подверженной в процессе нагружения возможности смещения как твердого тела, показаны преимущества созданного алгоритма по сравнению с формулировкой МКЭ в форме метода перемещений.

Физические и геометрические соотношения. Деформации и напряжения в произвольной точке оболочки связаны законом Гука

$\underset{3\times 1}{\left\{{\epsilon }^z\right\}}=\underset{3\times 3}{\left[C\right]}\underset{3\times 1}{\left\{\sigma \right\}},$ (1)

где ${\underset{1\times 3}{\left\{{\epsilon }^z\right\}}}^T=\left\{{\epsilon }_{11}^z{\epsilon }_{22}^z\text{\hspace{0.33em}2}{\epsilon }_{12}^z\right\};$ ${\underset{3\times 1}{\left\{\sigma \right\}}}^T=\left\{{\sigma }_{11}{\sigma }_{22}{\sigma }_{33}\right\}.$

На основе гипотезы прямой нормали деформации и напряжения определяются через параметры срединной поверхности

$\underset{3\times 1}{\left\{{\epsilon }^z\right\}}=\underset{3\times 6}{\left[G_{\epsilon }\right]}\underset{6\times 1}{\left\{\epsilon \right\}};\underset{3\times 1}{\left\{\sigma \right\}}=\underset{3\times 6}{\left[G_{\sigma }\right]}\underset{6\times 1}{\left\{N\right\}},$ (2)

где ${\underset{1\times 6}{\left\{\epsilon \right\}}}^T=\left\{{\epsilon }_{11}{\epsilon }_{22}\text{\hspace{0.33em}2}{\epsilon }_{12}{\aleph }_{11}{\aleph }_{22}\text{\hspace{0.33em}2}{\aleph }_{12}\right\};$ ${\underset{1\times 6}{\left\{N\right\}}}^T=\left\{N_{11}{N}_{22}{N}_{12}{M}_{11}{M}_{22}{M}_{12}\right\}.$

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

$v=e_{1}u+e_{2}v+e_{3}w;v_{,\alpha }=z_{\alpha }^k{e}_k;v_{,\mathrm{\alpha \beta }}=z_{\mathrm{\alpha \beta }}^k{e}_k.$ (3)

Деформации и искривления срединной поверхности определяются соотношениями Коши [8] в матричном виде

$\underset{6\times 1}{\left\{\epsilon \right\}}=\underset{6\times 3}{\left[L\right]}\underset{3\times 1}{\left\{\gamma \right\}},$ (4)

где ${\underset{1\times 3}{\left\{\gamma \right\}}}^T=\left\{uvw\right\}$ – строка перемещений точки срединной поверхности; $\underset{6\times 3}{\left[L\right]}$ – матрица алгебраических и дифференциальных операторов.

Элемент дискретизации. В соответствии с конечно-элементной моделью срединная поверхность оболочки представлялась совокупностью четырехугольных КЭ, узлы которых  $i,j,k,l$ располагались в вершинах четырехугольников. Искомыми узловыми неизвестными были выбраны продольные усилия и изгибающие моменты, а также компоненты вектора перемещения и их частные производные первого порядка.

Усилия и моменты точки срединной поверхности аппроксимируются через узловые значения с использованием билинейных функций:

$\underset{6\times 1}{\left\{N\right\}}=\underset{6\times 24}{\left[A\right]}\underset{24\times 1}{\left\{N_y\right\}},$ (5)

где ${\underset{1\times 24}{\left\{N_y\right\}}}^T=\left\{N_{11}^i{N}_{11}^j{N}_{11}^k{N}_{11}^l\mathrm{...}{N}_{12}^i{N}_{12}^j{N}_{12}^k{N}_{12}^l{M}_{11}^i\mathrm{...}{M}_{12}^l\right\}.$

Для интерполяции перемещений и их производных используются функции формы, элементами которых являются полиномы Эрмита третьей степени:

$\begin{array}{c}q={\underset{1\times 12}{\left\{\psi \right\}}}^T\underset{12\times 1}{\left\{q_y^L\right\}},\\ q_{,\alpha }={\underset{1\times 12}{\left\{{\psi }_{,\alpha }\right\}}}^T\underset{12\times 1}{\left\{q_y^L\right\},}\\ q_{,\alpha \beta }={\underset{1\times 12}{\left\{{\psi }_{,\alpha \beta }\right\}}}^T\underset{12\times 1}{\left\{q_y^L\right\}},\end{array}$ (6)

где ${\underset{1\times 12}{\left\{q_y^L\right\}}}^T=\left\{q^i{q}^j{q}^k{q}^l{q}_{,\xi }^i{q}_{,\xi }^j{q}_{,\xi }^k{q}_{,\xi }^l{q}_{,\eta }^i{q}_{,\eta }^j{q}_{,\eta }^k{q}_{,\eta }^l\right\}$ – строка узловых неизвестных в локальной системе координат, под символом q понимаются перемещения u, v, w.

Узловые неизвестные в локальной и глобальной системах координат связаны зависимостью

$\underset{12\times 1}{\left\{q_y^L\right\}}=\underset{12\times 12}{\left[D\right]}\underset{12\times 1}{\left\{q_y^G\right\}},$ (7)

где ${\underset{1\times 12}{\left\{q_y^G\right\}}}^T=\left\{q^i{q}^j{q}^k{q}^l{q}_{,S}^i{q}_{,S}^j{q}_{,S}^k{q}_{,S}^l{q}_{,\phi }^i{q}_{,\phi }^j{q}_{,\phi }^k{q}_{,\phi }^l\right\}.$

С использованием (6), (7) зависимость (4) приводится к виду

$\underset{6\times 1}{\left\{\epsilon \right\}}=\underset{6\times 36}{\left[B_1\right]}\underset{36\times 1}{\left\{U_y^G\right\}},$ (8)

где ${\underset{1\times 36}{\left\{U_y^G\right\}}}^T=\left\{{\underset{1\times 12}{\left\{u_y^G\right\}}}^T{\underset{1\times 12}{\left\{v_y^G\right\}}}^T{\underset{1\times 12}{\left\{w_y^G\right\}}}^T\right\}.$

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

Альтернативой аппроксимирующим выражениям (6) является векторная форма интерполяционной процедуры [7], которая предполагает, что в качестве объекта интерполяции выступает не какая-либо из компонент вектора перемещения, а непосредственно сам вектор в выбранной системе локальных ξ, η и глобальных S, φ координат:

${\underset{1\times 12}{\left\{v_y^L\right\}}}^T=\left\{v^i{v}^j{v}^k{v}^l{v}_{,\xi }^i{v}_{,\xi }^j{v}_{,\xi }^k{v}_{,\xi }^l{v}_{,\eta }^i{v}_{,\eta }^j{v}_{,\eta }^k{v}_{,\eta }^l\right\};$ (9)

${\underset{1\times 12}{\left\{v_y^G\right\}}}^T=\left\{v^i{v}^j{v}^k{v}^l{v}_{,S}^i{v}_{,S}^j{v}_{,S}^k{v}_{,S}^l{v}_{,\phi }^i{v}_{,\phi }^j{v}_{,\phi }^k{v}_{,\phi }^l\right\}.$ (10)

Аппроксимирующие выражения записываются для векторных величин:

$\begin{array}{c}v={\underset{1\times 12}{\left\{\psi \right\}}}^T\underset{12\times 1}{\left\{U_y^L\right\}}={\underset{1\times 12}{\left\{\psi \right\}}}^T\underset{12\times 12}{\left[D\right]}\underset{12\times 1}{\left\{U_y^L\right\}}=\\ ={\underset{1\times 12}{\left\{\pi \right\}}}^T\underset{12\times 36}{\left[EL\right]}\underset{36\times 36}{\left[PR\right]}\underset{36\times 1}{\left\{U_y^G\right\},}\\ v_{,\alpha }={\underset{1\times 12}{\left\{\pi \right\}}}^T\underset{12\times 36}{\left[EL\right]}\underset{36\times 36}{\left[PR\right]}\underset{36\times 1}{\left\{U_y^G\right\}},\\ v_{,\alpha \beta }={\underset{1\times 12}{\left\{\pi \right\}}}^T\underset{12\times 36}{\left[EL\right]}\underset{36\times 36}{\left[PR\right]}\underset{36\times 1}{\left\{U_y^G\right\}},\end{array}$ (11)

где $\underset{12\times 36}{\left[EL\right]}$ – матрица, элементами которой являются строки базисных векторов узловых точек.

С использованием выражений базисных векторов узловых точек через базисные векторы внутренней точки срединной поверхности конечного элемента выполняются координатные преобразования. В результате получаются аппроксимирующие соотношения для компонент векторов $v,v_{,\pm }$ и $v_{,{\pm }^2}$, с использованием которых зависимость (4) представляется матричным соотношением

$\underset{6\times 1}{\left\{\epsilon \right\}}=\underset{6\times 36}{\left[B_2\right]}\underset{36\times 1}{\left\{U_y^G\right\}}.$ (12)

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

$\begin{array}{c}{F}_S=\underset{V}{\int }{\underset{1\times 3}{\left\{{\sigma }_{\alpha \beta }\right\}}}^T\underset{3\times 1}{\left\{{\epsilon }^z\right\}}dV-\\ 0.5\underset{V}{\int }{\underset{1\times 3}{\left\{{\sigma }_{\alpha \beta }\right\}}}^T\underset{3\times 3}{\left[C_e\right]}\underset{3\times 1}{\left\{{\sigma }_{\alpha \beta }\right\}}dV-\\ 0.5\underset{S}{\int }{\underset{1\times 3}{\left\{\gamma \right\}}}^T\underset{3\times 1}{\left\{R\right\}}dS,\end{array}$ (13)

где $\left\{R\right\}$ – столбец заданной нагрузки.

При учете (2), (5), (8) и (12) функционал (13) можно представить в виде

$\begin{array}{c}{F}_S={\underset{1\times 24}{\left\{N_y\right\}}}^T\underset{V}{\int }{\underset{24\times 6}{\left[A\right]}}^T{\underset{6\times 3}{\left[G_{\sigma }\right]}}^T\underset{3\times 6}{\left[G_{\epsilon }\right]}\underset{6\times 36}{\left[B_{\alpha }\right]}dV\underset{36\times 1}{\left\{U_y^G\right\}}-\\ -0.5{\underset{1\times 24}{\left\{N_y\right\}}}^T\underset{V}{\int }{\underset{24\times 6}{\left[A\right]}}^T{\underset{6\times 3}{\left[G_{\sigma }\right]}}^T\underset{3\times 3}{\left[C_e\right]}\underset{3\times 6}{\left[G_{\sigma }\right]}\underset{6\times 24}{\left[A\right]}dV\underset{24\times 1}{\left\{N_y\right\}}-0.5{\underset{1\times 36}{\left\{U_y^G\right\}}}^T\underset{S}{\int }{\underset{36\times 3}{\left[GY\right]}}^T\underset{3\times 1}{\left\{R\right\}}dS,\end{array}$ (14)

где использовано выражение $\underset{3\times 1}{\left\{\gamma \right\}}={\underset{3\times 36}{\left[GY\right]}}^T\underset{36\times 1}{\left\{U_y^G\right\}}.$

Минимизация функционала (14) приводит к системе матричных уравнений

$\left\{\begin{array}{l}{\partial F_S/\partial \left\{N_y\right\}}^T\equiv -\underset{24\times 24}{\left[J\right]}\underset{24\times 1}{\left\{N_y\right\}}-\underset{24\times 36}{\left[H\right]}\underset{36\times 1}{\left\{U_y^G\right\}}=\mathrm{0,}\\ {\partial F_S/\partial \left\{U_y^G\right\}}^T\equiv {\underset{36\times 24}{\left[H\right]}}^T\underset{24\times 1}{\left\{N_y\right\}}-\underset{36\times 1}{\left\{f\right\}}=0.\end{array}\right.$ (15)

Из системы уравнений (15) можно выразить столбец силовых узловых неизвестных

$\underset{24\times 1}{\left\{N_y\right\}}={\underset{24\times 24}{\left[J\right]}}^{-1}\underset{24\times 36}{\left[H\right]}\underset{36\times 1}{\left\{U_y^G\right\}}$ (16)

и подставить его во второе уравнение (15):

${\underset{36\times 24}{\left[H\right]}}^T{\underset{24\times 24}{\left[J\right]}}^{-1}\underset{24\times 36}{\left[H\right]}\underset{36\times 1}{\left\{U_y^G\right\}}=\underset{36\times 1}{\left\{f\right\}},$ (17)

представляя (17) в виде

$\underset{36\times 36}{\left[K\right]}\underset{36\times 1}{\left\{U_y^G\right\}}=\underset{36\times 1}{\left\{f\right\}}.$ (18)

Матрица $\left[K\right]$ используется для компоновки глобальной матрицы податливости рассчитываемой оболочки вращения.

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

$r\left(x\right)=a+b\text{cos (}x/C).$ (19)

Входящие в (19) геометрические параметры имели следующие значения: a = 0.9 м; b = 0.4 м; c = 0.24 м. Осевая координата x изменялась в пределах 0 ≤ x ≤ 0.24π м. На левом торце оболочка первоначально имела шарнирное опирание, правый торец был свободен от закреплений. Вследствие наличия осевой симметрии оболочка моделировалась лентой элементов дискретизации, ориентированной вдоль меридиана оболочки (рис. 1).

 

Рис. 1. Расчетная схема оболочки вращения.

 

Прочие исходные данные были приняты равными:

E = 2 · 10⁵ МПа; ν = 0.3; q = 0.2 МПа; t = 0.01 м. В процессе выполнения численного эксперимента конечно-элементные расчеты осуществлены как в смешанной формулировке, так и в формулировке метода перемещений.

Результаты вычислений на основе скалярного варианта аппроксимации МКЭ представлены в табл. 1, где приведены значения нормальных напряжений в сечениях торцов оболочки при различных сетках дискретизации.

 

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

Координаты точек x, м; θ, рад.	Напряжения, МПа	Сетка узлов						Решение из условия равновесия
		МКЭ в форме метода перемещений			Смешанная формулировка МКЭ
		2 × 49	2 × 61	2 × 81	2 × 49	2 × 61	2 × 81
x = 0.0; θ = 0.0	${\mathrm{\sigma }}_{ss}^{in}$	11.57	11.38	11.23	11.02	11.04	11.05	–
	${\mathrm{\sigma }}_{ss}^{midl}$	11.25	11.16	11.11	11.068	11.071	11.074	11.077
	${\mathrm{\sigma }}_{ss}^{out}$	10.92	10.94	10.98	11.11	11.10	11.09	–
	${\mathrm{\sigma }}_{\theta \theta }^{in}$	−81.20	−81.27	−81.31	−81.35	−81.36	−81.36	–
	${\mathrm{\sigma }}_{\theta \theta }^{midl}$	−81.30	−81.33	−81.35	−81.34	−81.35	−81.35	–
	${\mathrm{\sigma }}_{\theta \theta }^{out}$	−81.40	−81.40	−81.39	−81.32	−81.34	−81.35	–
x = 0.24π; θ = 0.0	${\mathrm{\sigma }}_{ss}^{in}$	−0.291	−0.16	−0.074	−0.031	−0.021	−0.0123	–
	${\mathrm{\sigma }}_{ss}^{midl}$	−0.136	−0.073	−0.032	−0.0012	−0.0007	−0.0004	0.00
	${\mathrm{\sigma }}_{ss}^{out}$	0.021	0.015	0.0098	0.029	0.019	0.0115	–
	${\mathrm{\sigma }}_{\theta \theta }^{in}$	8.32	8.37	8.40	8.39	8.40	8.41	–
	${\mathrm{\sigma }}_{\theta \theta }^{midl}$	8.37	8.40	8.41	8.40	8.41	8.42	–
	${\mathrm{\sigma }}_{\theta \theta }^{out}$	8.42	8.42	8.43	8.41	8.41	8.42	–

 

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

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

${\sigma }_{ss}^{midl}=\frac{\pi {\left(a+b\right)}^2-\pi {\left(a-b\right)}^2}{2\pi \left(a+b\right)}\cdot \frac{q}{t}=\frac{{1.3}^2-{0.5}^2}{2\cdot 1.3}\cdot \frac{0.2}{0.01}=11.077$МПа.

Кроме того, очевидно, что на правом незагруженном торце меридиональные напряжения должны быть равны нулю (${\mathrm{\sigma }}_{ss}^{midl}=0.00$МПа).

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

Если шарнирные опоры на левом торце заменить на упругоподатливые, оболочка получит возможность смещаться вдоль осевой координаты х как абсолютно твердое тело. Величина данного жесткого смещения обратно пропорциональна жесткости пружинных опор. Очевидно, что наличие жестких смещений не должно оказывать никакого влияния на параметры НДС рассчитываемой оболочки. Для того, чтобы оценить влияние жестких смещений на величины контролируемых прочностных параметров НДС оболочки, были выполнены расчеты МКЭ при использовании приведенных вариантов интерполяционных процедур компонент вектора перемещения. В табл. 2 приведены результаты численного эксперимента по расчету НДС оболочки, допускающей смещение как твердого тела, полученные при реализации МКЭ в формулировке метода перемещений. Использование скалярной интерполяции компонент вектора перемещения (табл. 2) приводит к существенным погрешностям вычислений напряжений даже при незначительных смещениях оболочки как абсолютно твердого тела. При более значительных величинах жестких смещений оболочки погрешность вычислений прочностных параметров НДС резко возрастает до совершенно неприемлемых значений. Так, начиная с величины жесткого смещения вдоль оси Ox, равной 0.18 м и более, значения σ_ss на концевом сечении возрастают на несколько порядков, а σ_θθ вообще изменяют свой знак.

 

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

Координаты точек x, м; θ, рад.	Напряжения, МПа	Способы интерполяции кинематических искомых величин										Решение из условия равновесия
		Скалярный способ					Векторный способ
		Величина смещения оболочки как твердого тела, м
		0.00	0.09	0.18	0.377	0.754	0.00	0.09	0.18	0.377	0.754
x = 0.0; θ = 0.0	${\mathrm{\sigma }}_{ss}^{in}$	11.38	14.34	17.31	23.73	36.09	9.49	9.49	9.49	9.49	9.49	–
	${\mathrm{\sigma }}_{ss}^{midl}$	11.16	25.29	39.42	70.02	128.9	11.59	11.59	11.59	11.59	11.59	11.077
	${\mathrm{\sigma }}_{ss}^{out}$	10.94	36.24	61.53	116.3	221.6	13.69	13.69	13.69	13.69	13.69	–
	${\mathrm{\sigma }}_{\theta \theta }^{in}$	−81.27	−82.25	−83.24	−85.38	−89.49	−80.0	−80.0	−80.0	−80.0	−80.0	–
	${\mathrm{\sigma }}_{\theta \theta }^{midl}$	−81.33	−78.97	−76.61	−71.49	−61.66	−79.41	−79.41	−79.41	−79.41	−79.41	–
	${\mathrm{\sigma }}_{\theta \theta }^{out}$	−81.40	−75.69	−69.98	−57.60	−33.82	−78.80	−78.80	−78.80	−78.80	−78.80	–
x = 0.24π; θ = 0.0	${\mathrm{\sigma }}_{ss}^{in}$	−0.16	−25.34	−50.52	−105.1	−209.9	−0.108	−0.108	−0.108	−0.108	−0.108	–
	${\mathrm{\sigma }}_{ss}^{midl}$	−0.073	−14.09	−28.11	−58.48	−116.9	0.033	0.033	0.033	0.033	0.033	0.00
	${\mathrm{\sigma }}_{ss}^{out}$	−0.015	−2.84	−5.7	−11.89	−23.79	0.174	0.174	0.174	0.174	0.174	–
	${\mathrm{\sigma }}_{\theta \theta }^{in}$	8.37	−1.85	−12.07	−34.22	−76.78	8.51	8.51	8.51	8.51	8.51	–
	${\mathrm{\sigma }}_{\theta \theta }^{midl}$	8.40	1.52	−5.35	−20.24	−48.86	8.55	8.55	8.55	8.55	8.55	–
	${\mathrm{\sigma }}_{\theta \theta }^{out}$	8.42	4.90	1.37	−6.26	−20.94	8.59	8.59	8.59	8.59	8.59	–

 

Для корректного анализа НДС оболочек вращения в криволинейной системе координат следует использовать векторную аппроксимацию компонент вектора перемещения (9)–(12). Такой вид интерполяционной процедуры позволил корректным образом в полной мере учесть жесткие смещения оболочки. Анализируя соответствующие колонки табл. 2, можно констатировать, что при любой величине жестких смещений численные значения напряжений при реализации векторной формы интерполяции остаются абсолютно стабильными, что является необходимым условием корректности алгоритма анализа НДС оболочек вращения, допускающих жесткие смещения.

С целью визуализации результатов (табл. 2) на рис. 2 представлены графики изменения нормальных напряжений на срединной поверхности в зависимости от величины жесткого смещения оболочки вдоль оси Ox. Кривыми 1–4 показаны изменения значений напряжений σ_ss и σ_θθ при общепринятой в МКЭ форме интерполяционной процедуры (6) (кривые 1 и 2 отражают изменения напряжений в опорном сечении, а кривые 3 и 4 – в концевом сечении), а кривыми 5–8 – при использовании векторной формы интерполяции компонент вектора перемещений и их производных (9)–(12) (кривые 5 и 6 относятся к опорному сечению, а кривые 7 и 8 – к концевому сечению).

 

Рис. 2. Значения напряжений в формулировке метода перемещений с использованием скалярного (линии 1–4) и векторного (линии 5–8) способов интерполяции; линии 1, 3, 5, 7 – σ_ss, линии 2, 4, 6, 8 – σ_θθ.

 

Кривые 1–4 (рис. 2) имеют весьма существенный градиент, который нарастает по мере увеличения жесткого смещения, что является недопустимым. А кривые 5–8 представляют собой строго горизонтальные линии, что указывает на то, что значения напряжений не меняются даже при значительных величинах жесткого смещения оболочки при использовании векторной формы интерполяции компонент вектора перемещения.

Результаты расчета (табл. 3) показывают, что при использовании скалярной формы интерполяции значения напряжений σ_ss в опорном сечении практически не меняются, а значения σ_θθ возрастают по мере увеличения жесткого смещения, но отмеченная погрешность составляет 16%. В концевом сечении меридиональные напряжения σ_ss, несмотря на жесткие смещения, демонстрируют стабильность и остаются весьма близкими к нулевым значениям, что соответствует статическому условию равновесия незагруженного торца. Кольцевые же напряжения σ_θθ на правом конце существенно снижаются (до 20 раз) по мере увеличения жестких смещений, а при величине жесткого смещения 0.754 м меняют свой знак.

 

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

Координаты точек x, м; θ, рад.	Напряжения, МПа	Способы интерполяции кинематических искомых величин										Решение из условия равновесия
		Скалярный способ					Векторный способ
		Величина смещения оболочки как твердого тела, м
		0.00	0.09	0.18	0.377	0.754	0.00	0.09	0.18	0.377	0.754
x = 0.0; θ = 0.0	${\mathrm{\sigma }}_{ss}^{in}$	11.07	11.03	11.02	11.00	10.96	11.04	11.04	11.04	11.04	11.04	–
	${\mathrm{\sigma }}_{ss}^{midl}$	11.07	11.07	11.07	11.07	11.07	11.07	11.07	11.07	11.07	11.07	11.077
	${\mathrm{\sigma }}_{ss}^{out}$	11.10	11.11	11.12	11.14	11.18	11.11	11.11	11.11	11.11	11.11	–
	${\mathrm{\sigma }}_{\theta \theta }^{in}$	−81.36	−82.94	−84.53	−87.97	−94.58	−81.30	−81.30	−81.30	−81.30	−81.30	–
	${\mathrm{\sigma }}_{\theta \theta }^{midl}$	−81.35	−82.93	−84.52	−87.95	−94.55	−81.29	−81.29	−81.29	−81.29	−81.29	–
	${\mathrm{\sigma }}_{\theta \theta }^{out}$	−81.34	−82.92	−84.50	−87.93	−94.52	−81.28	−81.28	−81.28	−81.28	−81.28	–
x = 0.24π; θ = 0.0	${\mathrm{\sigma }}_{ss}^{in}$	−0.021	−0.053	−0.085	−0.154	−0.287	−0.021	−0.021	−0.021	−0.021	−0.021	–
	${\mathrm{\sigma }}_{ss}^{midl}$	−0.001	0.000	0.001	0.002	0.004	−0.001	−0.001	−0.001	−0.001	−0.001	0.00
	${\mathrm{\sigma }}_{ss}^{out}$	0.019	0.053	0.086	−0.158	0.296	0.019	0.019	0.019	0.019	0.019	–
	${\mathrm{\sigma }}_{\theta \theta }^{in}$	8.40	6.49	4.58	0.443	7.515	8.41	8.41	8.41	8.41	8.41	–
	${\mathrm{\sigma }}_{\theta \theta }^{midl}$	8.41	6.51	4.61	0.495	−7.418	8.42	8.42	8.42	8.42	8.42	–
	${\mathrm{\sigma }}_{\theta \theta }^{out}$	8.42	6.53	4.64	0.546	−7.321	8.42	8.42	8.42	8.42	8.42	–

 

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

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

 

Рис. 3. Значения напряжений в формулировке смешанного метода с использованием скалярного (линии 1–4) и векторного (линии 5–8) способов интерполяции; линии 1, 3, 5, 7 – σ_ss, линии 2, 4, 6, 8 – σ_θθ.

 

Кривые 1 и 3 (рис. 3), которые отражают изменение σ_ss при общепринятой в МКЭ форме интерполяции, – горизонтальные, а кривые 2 и 4, отражающие изменение σ_θθ при этой же форме интерполяции, имеют градиент, нарастающий при увеличении величины жесткого смещения. При использовании векторной формы интерполяционной процедуры (кривые 5–8) наблюдается абсолютная стабильность напряжений (кривые 5–8 строго горизонтальны) при всех рассматриваемых значениях величины смещения оболочки как жесткого целого, что подтверждает преимущества такой формы интерполяционной процедуры при анализе НДС оболочек вращения, в том числе допускающих жесткие смещения.

Пример расчета 2. С целью адекватной оценки точности вычисления прочностных параметров НДС оболочек вращения при использовании разработанного алгоритма на основе смешанного МКЭ с векторной аппроксимацией перемещений был выполнен сравнительный анализ конечно-элементных решений, полученных с применением разработанного алгоритма и вычислительного комплекса ANSYS. Первоначально параметр в формуле (33) был принят равным 0.10 м (0 ≤ x ≤ 0.10π м), затем параметру последовательно присваивались значения c = 0.08 м; c = 0.06 м; c = 0.05 м. Была использована расчетная схема (рис. 1). Прочие геометрические и физические значения исходных данных соответствовали приведенным выше величинам. Результаты расчетов представлены в табличной форме (табл. 4–7).

 

Таблица 4. Значения нормальных напряжений при использовании вычислительного комплекса ANSYS и в смешанной формулировке при c = 0.10 м

Координаты точек x, м; θ, рад.	Напряжения, МПа	Варианты расчета									Решение из условия статического равновесия
		ANSYS Shell 181					Смешанный МКЭ с тензорно-векторной интерполяцией
		Сетка узлов
		21 × 31	21 × 61	21 × 101	21 × 151	21 × 201	2 × 31	2 × 61	2 × 81	2 × 101
x = 0.0; θ = 0.0	${\mathrm{\sigma }}_{ss}^{in}$	16.53	13.07	11.92	11.58	11.63	11.62	11.21	11.15	11.12	–
	${\mathrm{\sigma }}_{ss}^{midl}$	10.07	10.88	11.125	11.205	11.42	10.98	11.05	11.06	11.07	11.077
	${\mathrm{\sigma }}_{ss}^{out}$	3.61	8.69	10.33	10.83	11.21	10.35	10.89	10.97	11.01	–
	${\mathrm{\sigma }}_{\theta \theta }^{in}$	−349.22	−379.92	−394.3	−401.88	−408.7	−402.00	−403.3	−403.5	−403.5	–
	${\mathrm{\sigma }}_{\theta \theta }^{midl}$	−356.65	−383.2	−395.7	−402.25	−407.7	−402.1	−403.3	−403.5	−403.6	–
	${\mathrm{\sigma }}_{\theta \theta }^{out}$	−364.07	−386.53	−397.01	−402.62	−406.7	−402.1	−403.4	−403.5	−403.6	–
x = 0.10π; θ = 0.0	${\mathrm{\sigma }}_{ss}^{in}$	−0.371	−0.216	−0.098	−0.048	−0.003	−0.109	−0.0310	−0.0179	−0.0117	–
	${\mathrm{\sigma }}_{ss}^{midl}$	0.325	0.157	0.075	0.034	0.031	−0.0045	−0.0011	−0.0006	−0.0004	0.00
	${\mathrm{\sigma }}_{ss}^{out}$	1.02	0.529	0.247	0.116	0.064	0.0995	0.0288	0.0167	0.0109	–
	${\mathrm{\sigma }}_{\theta \theta }^{in}$	1.8	1.79	−0.399	−1.61	−2.49	−3.131	−3.057	−3.042	−3.037	–
	${\mathrm{\sigma }}_{\theta \theta }^{midl}$	3.94	2.06	−0.402	−1.80	−2.95	−3.327	−3.103	−3.070	−3.055	–
	${\mathrm{\sigma }}_{\theta \theta }^{out}$	6.08	2.32	−0.405	−1.98	−3.41	−3.524	−3.152	−3.097	−3.072	–

 

Как следует из анализа значений напряжений (табл. 4) (c = 0.10 м) при существенном сгущении сетки КЭ, величины нормальных напряжений, полученные при использовании вычислительного комплекса ANSYS в опорном сечении, оказались примерно равными значениям напряжений, вычисленным при использовании разработанного алгоритма.

В концевом же сечении значения кольцевых напряжений, вычисленные в комплексе ANSYS при относительно редкой сетке КЭ, имеют противоположный знак значениям, полученным при достаточно густой сетке элементов дискретизации. Для достижения уровня точности вычисления кольцевых напряжений, соответствующего уровню точности разработанного алгоритма смешанного МКЭ с векторной аппроксимацией перемещений, вычислительному комплексу ANSYS понадобилось как минимум в 3–4 раза большее число КЭ в меридиональном направлении. Анализируя значения нормальных напряжений, вычисленных при использовании разработанного алгоритма, можно констатировать гораздо лучшую сходимость вычислительного процесса по сравнению со сходимостью, продемонстрированной комплексом ANSYS.

При уменьшении значения параметра до c = 0.08 м кривизна меридиана оболочки пропорционально возрастет. Результаты расчета при c = 0.08 м приведены в табл. 5.

 

Таблица 5. Значения нормальных напряжений при использовании вычислительного комплекса ANSYS и в смешанной формулировке при c = 0.08 м

Координаты точек x, м; θ, рад.	Напряжения, МПа	Варианты расчета									Решение из условия статического равновесия
		ANSYS Shell 181					Смешанный МКЭ с тензорно-векторной интерполяцией
		Сетка узлов
		21 × 31	21 × 61	21 × 101	21 × 151	21 × 201	2 × 31	2 × 61	2 × 81	2 × 121
x = 0.0; θ = 0.0	${\mathrm{\sigma }}_{ss}^{in}$	17.94	14.55	12.89	12.26	12.06	11.86	11.28	11.19	11.13	–
	${\mathrm{\sigma }}_{ss}^{midl}$	9.177	10.55	11.135	11.365	11.47	10.94	11.04	11.056	11.07	11.077
	${\mathrm{\sigma }}_{ss}^{out}$	0.414	6.55	9.38	10.47	10.88	10.01	10.80	10.92	11.01	–
	${\mathrm{\sigma }}_{\theta \theta }^{in}$	−449.1	−492.2	−514.05	−526.1	−532.4	−520.7	−522.1	−522.3	−522.4	–
	${\mathrm{\sigma }}_{\theta \theta }^{midl}$	−463.95	−500.45	−518.28	−528.1	−533.15	−520.7	−522.1	−522.3	−522.45	–
	${\mathrm{\sigma }}_{\theta \theta }^{out}$	−478.8	−508.7	−522.5	−530.1	−533.9	−520.7	−522.2	−522.4	−522.5	–
x = 0.08π; θ = 0.0	${\mathrm{\sigma }}_{ss}^{in}$	−0.733	−0.415	−0.21	−0.12	−0.09	−0.125	−0.0346	−0.0199	−0.009	–
	${\mathrm{\sigma }}_{ss}^{midl}$	0.379	0.223	0.118	0.001	0.025	−0.0046	−0.0011	−0.0006	−0.0003	0.00
	${\mathrm{\sigma }}_{ss}^{out}$	1.49	0.86	0.445	0.123	0.14	0.115	0.0323	0.0186	0.0084	–
	${\mathrm{\sigma }}_{\theta \theta }^{in}$	3.62	−3.25	−6.91	−9.01	−10.1	−12.36	−12.47	−12.49	−12.50	–
	${\mathrm{\sigma }}_{\theta \theta }^{midl}$	5.29	−2.165	−6.37	−8.88	−10.25	−12.84	−12.59	−12.56	−12.53	–
	${\mathrm{\sigma }}_{\theta \theta }^{out}$	6.95	−1.08	−5.83	−8.75	−10.4	−13.33	−12.72	−12.63	−12.56	–

 

Анализируя табличные данные, можно отметить исключительно быструю сходимость вычислительного процесса разработанного алгоритма по сравнению с ANSYS. Так, значения меридиональных и кольцевых напряжений на срединной поверхности в опорном сечении при сетке 2 × 31 и 2 × 121 изменяются всего лишь на 1.3 и 0.33% соответственно. В концевом сечении данный процент по кольцевым напряжениям составляет 2.5. Меридиональные напряжения в концевом сечении оказались достаточно близки к нулевым значениям даже при редкой сетке КЭ. При использовании программного комплекса ANSYS процент расхождения при сетке узлов 21 × 61 и 21 × 201 (сетка 21 × 31 во внимание даже не принималась) в опорном сечении составил 8.0 и 6.1 для меридиональных и кольцевых напряжений на срединной поверхности соответственно. В концевом же сечении для кольцевых напряжений процент расхождения составил уже 78.9. Кроме того, даже при густой сетке КЭ (21 × 201) кольцевые напряжения в концевом сечении, вычисленные в программном комплексе ANSYS, оказались заниженными по абсолютной величине на 20% по сравнению с кольцевыми напряжениями, полученными при использовании разработанного алгоритма.

По данным табл. 6, которые были получены при значении параметра c = 0.06 м, по сравнению с данными табл. 5 сходимость вычислительного процесса программного комплекса ANSYS ухудшилась, а проценты расхождений между значениями напряжений при сетке узлов 21 × 41 и 21 × 201 увеличились до 31.5 и 10.4 для меридиональных и кольцевых напряжений в опорном сечении соответственно. Кольцевые значения в концевом сечении, вычисленные при сетке 21 × 201 в комплексе ANSYS, оказались заниженными по сравнению с разработанным алгоритмом уже на 31.8%.

 

Таблица 6. Значения нормальных напряжений при использовании вычислительного комплекса ANSYS и в смешанной формулировке при c = 0.06 м

Координаты точек x, м; θ, рад.	Напряжения, МПа	Варианты расчета									Решение из условия статического равновесия
		ANSYS Shell 181					Смешанный МКЭ с тензорно-векторной интерполяцией
		Сетка узлов
		21 × 41	21 × 81	21 × 121	21 × 151	21 × 201	2 × 41	2 × 61	2 × 81	2 × 151
x = 0.0; θ = 0.0	${\mathrm{\sigma }}_{ss}^{in}$	17.79	15.62	14.5	14.25	13.63	11.68	11.35	11.23	11.12	–
	${\mathrm{\sigma }}_{ss}^{midl}$	8.21	10.305	11.24	11.725	11.99	10.94	11.01	11.04	11.07	11.077
	${\mathrm{\sigma }}_{ss}^{out}$	−1.37	4.99	7.98	9.2	10.35	10.20	10.68	10.85	11.01	–
	${\mathrm{\sigma }}_{\theta \theta }^{in}$	−609.0	−660.2	−682.7	−692.9	−703.8	−681.8	−682.1	−682.3	−682.4	–
	${\mathrm{\sigma }}_{\theta \theta }^{midl}$	−636.75	−678.2	−695.15	−702.65	−710.45	−681.5	−682.1	−682.3	−682.4	–
	${\mathrm{\sigma }}_{\theta \theta }^{out}$	−664.5	−696.2	−707.6	−712.4	−717.1	−681.2	−682.1	−682.3	−682.4	–
x = 0.06π; θ = 0.0	${\mathrm{\sigma }}_{ss}^{in}$	−0.84	−0.48	−0.326	−0.27	−0.23	−0.074	−0.034	−0.0194	−0.0056	–
	${\mathrm{\sigma }}_{ss}^{midl}$	0.325	0.195	0.113	0.075	0.031	−0.003	−0.0013	−0.0007	−0.0002	0.00
	${\mathrm{\sigma }}_{ss}^{out}$	1.49	0.87	0.552	0.42	0.29	0.068	0.0315	0.018	0.0052	–
	${\mathrm{\sigma }}_{\theta \theta }^{in}$	−2.86	−12.61	−16.67	−18.44	−20.25	−28.76	−29.12	−29.25	−29.37	–
	${\mathrm{\sigma }}_{\theta \theta }^{midl}$	0.70	−10.21	−15.24	−17.6	−20.09	−29.59	−29.50	−29.46	−29.43	–
	${\mathrm{\sigma }}_{\theta \theta }^{out}$	4.26	−7.81	−13.81	−16.69	−19.96	−30.43	−29.87	−29.67	−29.48	–

 

Если значение параметра уменьшить до c = 0.05 м, то кривизна меридиана оболочки возрастет, и погрешность значений напряжений, полученных с помощью программного комплекса ANSYS, возрастет. Это подтверждается анализом результатов (табл. 7) при значении параметра c = 0.05 м. Если в опорном сечении сходимость вычислительного процесса еще наблюдается (хотя ${\mathrm{\sigma }}_{ss}^{midl}$ отличаются на 12.4% от σ_ss по условиям равновесия при сетке узлов 21 × 201, а ${\mathrm{\sigma }}_{\theta \theta }^{midl}$ монотонно возрастают от 775.25 до 832.4 МПа), то в концевом сечении сходимость по кольцевым напряжениям явно неудовлетворительная. Кроме того, значения σ_θθ на внутренней и наружной поверхностях оболочки имеют разный знак, чего не наблюдалось при предыдущих значениях параметра с, равных 0.10, 0.08, 0.06 м. Также следует отметить, что при значении параметра с = 0.05 м нарушается выявленная при других значениях параметра с закономерность, которая заключается в том, что при уменьшении параметра с напряжения σ_θθ в концевом сечении возрастают по абсолютной величине. Анализируя правую часть табл. 7, в которой представлены значения напряжений, вычисленных при использовании разработанного алгоритма, можно наблюдать абсолютную сходимость вычислительного процесса при любом варианте сетки дискретизации и практическое совпадение меридиональных напряжений ${\mathrm{\sigma }}_{ss}^{midl}$ с напряжениями, вычисленными из условия статического равновесия оболочки.

 

Таблица 7. Значения нормальных напряжений при использовании вычислительного комплекса ANSYS и в смешанной формулировке при c = 0.05 м

Координаты точек x, м; θ, рад.	Напряжения, МПа	Варианты расчета									Решение из условия статического равновесия
		ANSYS Shell 181					Смешанный МКЭ с тензорно-векторной интерполяцией
		Сетка узлов
		21 × 61	21 × 101	21 × 151	21 × 181	21 × 201	2 × 61	2 × 81	2 × 101	2 × 151
x = 0.0; θ = 0.0	${\mathrm{\sigma }}_{ss}^{in}$	16.33	15.79	15.21	14.98	14.87	11.37	11.24	11.18	11.12	–
	${\mathrm{\sigma }}_{ss}^{midl}$	8.62	10.49	11.78	12.23	12.45	10.99	11.02	11.04	11.06	11.077
	${\mathrm{\sigma }}_{ss}^{out}$	0.919	5.19	8.34	9.47	10.02	10.60	10.80	10.90	11.0	–
	${\mathrm{\sigma }}_{\theta \theta }^{in}$	−740.8	−780.3	−805.65	−815.4	−820.4	−782.9	−782.9	−782.9	−782.9	–
	${\mathrm{\sigma }}_{\theta \theta }^{midl}$	−775.25	−805.15	−822.63	−829.05	−832.4	−782.7	−782.9	−782.9	−782.9	–
	${\mathrm{\sigma }}_{\theta \theta }^{out}$	−809.7	−830.0	−839.6	−842.7	−844.4	−782.5	−782.8	−782.9	−782.95	–
x = 0.05π; θ = 0.0	${\mathrm{\sigma }}_{ss}^{in}$	−1.17	−0.69	−0.44	−0.407	−0.39	−0.026	−0.0146	−0.0094	−0.0041	–
	${\mathrm{\sigma }}_{ss}^{midl}$	0.13	0.085	0.005	0.017	0.015	−0.0015	−0.0008	−0.0005	−0.0002	0.00
	${\text{σ}}_{ss}^{out}$	1.43	0.86	0.45	0.44	0.42	0.0227	0.0143	0.0084	0.0037	–
	${\mathrm{\sigma }}_{\theta \theta }^{in}$	6.89	0.602	−6.34	−5.17	−4.78	−33.27	−33.55	−33.67	−33.80	–
	${\mathrm{\sigma }}_{\theta \theta }^{midl}$	14.06	6.63	−1.7	−1.35	−1.22	−33.99	−33.45	−33.93	−33.90	–
	${\mathrm{\sigma }}_{\theta \theta }^{out}$	21.23	12.65	2.94	2.48	2.34	−34.71	−34.35	−34.18	−34.01	–

 

Выводы.

1. МКЭ в смешанной формулировке позволяет непосредственно получать внутренние силовые факторы (усилия и моменты) в интересующих узловых точках без дополнительных вычислительных процедур, обязательных при использовании МКЭ в форме метода перемещений.
2. Вычислительный алгоритм МКЭ в смешанной формулировке, предназначенный для анализа НДС тонких оболочек вращения, менее сложный по сравнению с алгоритмом МКЭ в форме метода перемещений из-за отсутствия необходимости организации вычислительной процедуры по нахождению численных значений вторых производных перемещений вдоль нормалей к срединной поверхности и дальнейшего их использования в соотношениях Коши для определения деформаций и напряжений в узловых точках.
3. При наличии жестких смещений у рассчитываемой оболочки МКЭ в смешанной формулировке демонстрируется лучшая, по сравнению с МКЭ в форме метода перемещений, устойчивость вычислительного процесса при использовании общепринятой интерполяционной процедуры перемещений как скалярных полей.
4. На основании выполненного сравнительного анализа результатов конечно-элементных решений, полученных с использованием комплекса ANSYS и разработанного алгоритма, можно сделать вывод, что при расчете пологих оболочек вращения контролируемые параметры НДС, вычисленные с помощью ANSYS и разработанного алгоритма, оказываются достаточно близкими по своим значениям. Но по мере увеличения кривизны меридиана оболочки погрешность вычислений в программном комплексе ANSYS возрастает и может превысить предельно допустимый уровень. Значения напряжений, полученные при использовании разработанного алгоритма расчета на основе смешанного МКЭ с векторной аппроксимацией перемещений, оказываются корректными даже при относительно редкой сетке узлов дискретизации, что позволяет эффективно исследовать НДС оболочек вращения при минимальных затратах времени и ограниченных ресурсах вычислительной техники.

Конфликт интересов. Авторы заявляют об отсутствии конфликта интересов.

Авторлар туралы

Ю. Клочков

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

Хат алмасуға жауапты Автор.
Email: klotchkov@bk.ru
Ресей, Волгоград

В. Пшеничкина

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

Email: klotchkov@bk.ru
Ресей, Волгоград

А. Николаев

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

Email: klotchkov@bk.ru
Ресей, Волгоград

С. Марченко

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

Email: klotchkov@bk.ru
Ресей, Волгоград

О. Вахнина

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

Email: klotchkov@bk.ru
Ресей, Волгоград

М. Клочков

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

Email: klotchkov@bk.ru
Ресей, Волгоград

Әдебиет тізімі

Қосымша файлдар