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

Cover Page

Cite item

Full Text

Abstract

A three-field finite element of a quadrangular shape of a thin shell with nodal unknowns in the form of: displacements and their first derivatives has been developed; deformations and curvatures of the median surface; forces and moments of the middle surface.

The approximation of the required quantities was carried out in two versions. In the first version, the components of the displacement vector and the components of the strain and curvature tensors, as well as the force and moment tensors, were approximated using traditional shape functions as components of scalar fields. In the second version, tensor quantities were approximated through the corresponding tensors of nodal points, and only after coordinate transformations based on the relations of the used curvilinear coordinate system were approximating expressions for the components of the corresponding tensors obtained.

Specific examples show the effectiveness of using the second version of approximating expressions in shell calculations.

Full Text

  1. Введение

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

R0=xi+rx, θsinθj+rx, θcosθk, (2.1)

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

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

rx, θ=bca2x2ab2cos2θ+c2sin2θ, (2.2)

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

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

rθ=bcb2cos2θ+c2sin2θ (2.3)

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

a10=Rx0, a20=Rθ0, a0=a10×a20/a0, (2.4)

где a0=a10a10a20a20a10a202.

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

g10=Rx0ζ, g20=Rθ0ζ, (2.5)

где R0ζ=R0+ζa0.

В процессе деформирования оболочки под действием приложенной внешней нагрузки точки M0 и M займут новые положения M и Mζ, определяемые векторами перемещений v и V соответственно

v=v1a10+v2a20+va0, V=v+ζaa0 (2.6)

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

a=a1×a2/a, (2.7)

где a1=Rx=R0+vx, a2=R,θ=R0+v,θ, a=a1a1a2a2a1a22.

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

g1=Rxζ=R0ζ+Vx=g10+vx+ζaa0xg2=Rθζ=R0ζ+Vθ=g20+vθ+ζaa0θ (2.8)

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

gαβ0=gα0gβ0, gαβ=gαgβ, (2.9)

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

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

εαβζ=gαβgαβ0/2. (2.10)

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

εαβζ3×1=εαβ3×1+ζαβ3×1, (2.11)

где εαβζ1×3T=ε11ζε22ζ2ε12ζ, εαβ1×3T=ε11ε222ε12, αβ1×3T=1122212 – деформации и искривления в точках Mζ и M оболочки.

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

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

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

x=φ1×4Txy4×1, θ=φ1×4Tθy4×1, (3.1)

где xy1×4T=xi xj xk xl, θy1×4T=θi θj θk θl.

В качестве искомых кинематических и силовых параметров конечного элемента выбираются следующие столбцы узловых неизвестных: vyL1×36T=vy1L1×12Tvy2L1×12TvyL1×12T – строки узловых перемещений, каждая из которых имеет вид qyL1×12T=qi qj qk ql qξi...qξl qηi...qηl, где под q понимается компонента вектора перемещения v1, v2 или v; εy1×24T=εy1×12Ty1×12T – строки узловых деформаций и искривлений срединной поверхности λy1×12T=λ11i λ11j λ11k λ11l λ22i...λ22l 2λ12i...2λ12l, где под λαβ понимается деформация εαβ или искривления ℵαβ срединной поверхности; NMy1×24T=Ny1×12TMy1×12T – строки узловых продольных сил и изгибающих моментов Ty1×12T=T11i T11j T11k T11l T22i...T22l T12i...T12l. Здесь под Tαβ понимается контравариантная компонента тензора внутренних усилий или тензора моментов.

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

q=ψ1×12TqyL12×1, λαβ=φ1×4Tλαβy4×1, Tαβ=φ1×4TTyαβ4×1, (3.2)

где ψ1×12T – матрица-строка, содержащая произведения полиномов Эрмита третьей степени [24, 25].

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

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

ε6×1=H16×24εy24×1, NM6×1=H16×24NMy24×1, v3×1=A13×36vy36×1, εk6×1=B16×36vy36×1, (3.3)

где ε1×6T=εαβT1×3αβT1×3 – строка деформаций и искривлений внутренней точки конечного элемента; NM1×6T=N11 N22 N12 M11 M22 M12 – строка усилий и моментов внутренней точки конечного элемента; εk1×6 – строка деформаций и искривлений внутренней точки, определяемых функциями перемещений.

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

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

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

ε~=φ1×4Tε~y4×1, ~=φ1×4T~y4×1, N~=φ1×4TN~y4×1, M~=φ1×4TM~y4×1, v=ψ1×12TvyL12×1, (3.4)

где ε~y1×4T=ε~i ε~j ε~k ε~l; ~y1×4T=~i ~j ~k ~l – строки узловых тензоров деформаций и искривлений срединной поверхности; N~y1×4T=N~i N~j N~k N~lM~y1×4T=M~i M~j M~k M~l – строки узловых тензоров продольных сил и изгибающих моментов; vyL1×12T=vi vj vk vl vξi...vξl vηi...vηl – строка узловых векторов перемещения и их производных первого порядка.

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

ε~=ε11a01a01+ε22a02a02+2ε12a01a02=a0αβ1×3Tεαβ3×1, ~=a0αβ1×3Tαβ3×1N~=N11a10a10+N22a20a20+N12a10a20=aαβ01×3TNαβ3×1, M~=aαβ01×3TMαβ3×1ε~y4×1=D4×12εαβy12×1, ~y4×1=D4×12αβy12×1, N~y4×1=L4×12Nyαβ12×1, M~y4×1=L4×12Myαβ12×1, (3.5)

где

L4×12=a10a10ia20a20ia10a20i000000000000a10a10ja20a20ja10a20j000000000000a10a10ka20a20ka10a20k000000000000a10a10la20a20la10a20l

матрица D4×12 имеет вид L4×12, но ее компоненты, отличные от нуля, являются диадными произведениями контравариантных векторов базисов узловых точек;

εαβy1×12T=ε11i ε22i 2ε12i1×3ε11j ε22j 2ε12j1×3ε11k ε22k 2ε12k1×3ε11l ε22l 2ε12l1×3αβy1×12T=11i 22i 212i1×311j 22j 212j1×311k 22k 212k1×311l 22l 212l1×3Nyαβ1×12T=N11i N22i N12i1×3N11j N22j N12j1×3N11k N22k N12k1×3N11l N22l N12l1×3Myαβ1×12T=M11i M22i M12i1×3M11j M22j M12j1×3M11k M22k M12k1×3M11l M22l M12l1×3

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

aα03×1=m03×3i3×1, aα0ρ3×1=m0ρ3×3i3×1, (3.6)

где aα01×3T=a10 a10 a0, i1×3T=i j k, aα0ρ1×3T=a10ρ a20ρ a0ρ. Верхний индекс ρ обозначает конкретный из узлов КЭ i, j, k, l.

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

aα0ρ3×1=m0ρ3×3m03×3-1aα03×1=ωρ3×3aα03×1 (3.7)

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

aαβ0ρ3×1=b0ρ3×3aαβ03×1, a0αβρ3×1=d0ρ3×3a0αβ3×1, (3.8)

где aαβ0ρ1×3T=a10ρa10ρa20ρa20ρa10ρa20ρ, a0αβρ1×3T=a01a01a02a02a01a02.

Принимая во внимание (3.8), матрица L примет вид

L4×12=aαβ0Tb0iT000000000000aαβ0Tb0jT000000000000aαβ0Tb0kT000000000000aαβ0Tb0lT (3.9)

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

a0αβ1×3Tεαβ3×1=a0αβ1×3Tϕ1d0iT3×3ϕ2d0jT3×3ϕ3d0kT3×3ϕ4d0lT3×3εαβy12×1a0αβ1×3Tαβ3×1=a0αβ1×3Tϕ1d0iT3×3ϕ2d0jT3×3ϕ3d0kT3×3ϕ4d0lT3×3αβy12×1aαβ01×3TNαβ3×1=aαβ01×3Tϕ1b0iT3×3ϕ2b0jT3×3ϕ3b0kT3×3ϕ4b0lT3×3Nyαβ12×1aαβ01×3TMαβ3×1=aαβ01×3Tϕ1b0iT3×3ϕ2b0jT3×3ϕ3b0kT3×3ϕ4b0lT3×3Myαβ12×1 (3.10)

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

ε11ε222ε12=ϕ1d0i3×3Tϕ2d0j3×3Tϕ3d0k3×3Tϕ4d0l3×3Tεαβy12×11122212=ϕ1d0i3×3Tϕ2d0j3×3Tϕ3d0k3×3Tϕ4d0l3×3Tαβy12×1

N11N22N12=ϕ1b0i3×3Tϕ2b0j3×3Tϕ3b0k3×3Tϕ4b0l3×3TNyαβ12×1M11M22M12=ϕ1b0i3×3Tϕ2b0j3×3Tϕ3b0k3×3Tϕ4b0l3×3TMyαβ12×1 (3.11)

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

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

ε6×1=H26×24εy24×1, NM6×1=H26×24NMy24×1, v3×1=A23×36vy36×1, ξk6×1=B26×36vy36×1. (3.12)

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

При общепринятой форме интерполяционной процедуры (3.2) отдельная компонента тензора, например, компонента тензора искривлений ℵ11 является функцией узловых значений только этой же самой компоненты, т.е. 11=f11i 11j 11k 11l и не зависит от узловых значений остальных компонент данного тензора 22ρ, 12ρ.

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

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

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

L=FNM1×6Tε6×1dFFε1×6TNM6×1C6×6εk6×1dFFU1×3TP3×1dF, (4.1)

где NM1×6T=N11 N22 N12 M11 M22 M12, ε1×6T=ε11 ε22 2ε12 11 22 212 – матрицы-строки продольных сил и изгибающих моментов, а также деформационных параметров в точке срединной поверхности оболочки; {εℵk} – столбец, элементы которого определяются соотношениями Коши; U1×3T=v1 v2 vP1×3T=p1 p2 p3 – матрицы-строки компонент вектора перемещения и вектора внешней поверхностной нагрузки.

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

NM6×1=C6×6ε6×1 (4.2)

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

L=NMy1×24TFHγ24×6TC6×61Hγ6×24dFNMy24×1εy1×24TFHγ24×6THγ6×24dFNMy24×1++εy1×24TFHγ24×6TC6×6Bγ6×36dFPR36×36vyG36×1vyG1×36TPR36×36TFAγ36×3P3×1dF, (4.3)

где матрицы Hγ6×24, Bγ6×36 и Aγ36×3 компонуются на основе интерполяционных выражений (3.3) или (3.12), нижний индекс γ принимает значения 1, 2, причем γ = 1 соответствует использованию общепринятой покомпонентной интерполяционной процедуре (3.2), а γ = 2 соответствует применению разработанной тензорно-векторной форме интерполяции искомых величин (3.11) в трехпольном варианте МКЭ.

Входящая в (4.3) матрица PR36×36 определяет связь между столбцами кинематических узловых неизвестных в локальной и глобальной системах координат.

Осуществляя минимизацию (4.3) по кинематическим vyGT, {εℵy}T и силовым {NMy}T искомым неизвестным, можно записать систему матричных уравнений

LvyGT=PR36×36TFB36×6TC6×6THγ6×24dFεy24×1PR36×36TFAγ36×3P3×1dF=0LεyT=FHγ24×6TC6×6B6×36dFPR36×36vyG36×1FHγ24×6THγ6×24dFNMy24×1=0LNMyT=FHγ24×6TC6×61Hγ6×24dFNMy24×1FHγ24×6THγ6×24dFεy24×1=0 (4.4)

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

S36×24Tεy24×1f36×1=0S24×36vyG36×1G24×24NMy24×1=0W24×24NMy24×1G24×24Tεy24×1=0, (4.5)

где S24×36=FHγ24×6TC6×6B6×36dFPR36×36, G24×24=FHγ24×6THγ6×24dF, W24×24=FHγ24×6TC6×61Hγ6×24dF.

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

εy24×1=GT24×241W24×24NMy24×1, NMy24×1=G24×241S24×36vyG36×1εy24×1=GT24×241W24×24G24×241S24×36vyG36×1 (4.6)

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

K36×36vyG36×1=f36×1, (4.7)

где K36×36=S36×24TG124×24TW24×24G24×24-1S24×36 – искомая матрица жесткости размером 36 × 36 конечного элемента, используемого для построения дискретной модели оболочки.

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

Полученные в результате решения методом Гаусса глобальной СЛАУ столбцы vyG посредством (4.6), позволяют вычислять требуемые для анализа НДС оболочки столбцы продольных сил и изгибающих моментов {εℵy}, а также столбцы деформаций и искривлений {NMy} в любой интересующей исследователя узловой точке рассчитываемой конструкции.

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

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

 

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

 

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

Nф22=qL=0.11=0.100 Нσθθmidl=Nф22/Lh=0.100 Н/1 см0.1 см=1.000 Н/см2Mф22=0.000 Нсм

 

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

Координаты точек, y, м; θ, рад.

σ, 10–2 МПа; M, Н · см; N, Н

Способ интерполяции

Значения M22, N22, σθθmidl из условия равновесия

Общепринятый

Тензорно-векторный

Сетка узлов

121 × 2151 × 2181 × 2211 × 2121 × 2151 × 2181 × 2211 × 2

y = 0.0; θ = 0.0

σθθin

971.4

971.4

971.4

971.4

971.3

971.3

971.4

971.4

σθθout

–977.7

–977.7

–977.7

–977.75

–977.7

–977.7

–977.7

–977.7

M22

–1.624

–1.624

–1.624

–1.624

–1.624

–1.624

–1.624

–1.624

σxxin

291.4

291.4

291.4

291.4

291.4

291.4

291.4

291.4

σxxout

–293.3

–293.3

–293.3

–293.3

–293.3

–293.3

–293.3

–293.3

y = 0.5; θ=π2

σθθin

0.319

–0.405

–0.703

–0.840

0.949

–0.091

–0.524

–0.728

σθθout

–2.151

–1.508

–1.247

–1.129

–2.834

–1.849

–1.441

–1.250

σθθmidl

–0.974

–0.983

–0.988

–0.992

–1.035

–1.013

–1.005

–1.002

–1.000

M22

–0.002

–0.001

–0.001

–0.0003

–0.003

–0.0016

–0.001

–0.0005

0.000

N22

–0.097

–0.098

–0.099

–0.099

–0.1035

–0.101

–0.100

–0.100

–0.100

σxxin

1.825

0.612

0.157

–0.041

2.135

0.993

0.475

0.214

σxxout

0.540

0.061

–0.110

–0.183

0.664

0.258

0.067

–0.039

 

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

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

 

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

 

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

 

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

Координаты точек, y, м; θ, рад.

σ, 10–2 МПа; M, Н · см; N, Н

Величина жесткого смещения, м

Значения M22, N22, σθθmidl из условия равновесия

0.00

0.05

0.10

0.25

0.50

0.75

Общепринятый способ интерполяции

y = 0.0; θ = 0.0σθθin

971.4

971.2

971.1

970.6

969.8

969.0

σθθout

–977.75

–977.6

–977.4

–976.9

–976.1

–975.3

M22

–1.624

–1.624

–1.624

–1.623

–1.622

–1.620

σxxin

291.4

291.4

291.35

291.25

291.1

290.9

σxxout

–293.3

–293.3

–293.2

–293.0

–292.7

–292.4

y = 0.5; θ=π2σθθin

–0.840

–0.246

0.346

2.107

4.994

7.825

σθθout

–1.129

–0.770

–0.414

0.648

2.390

4.098

σθθmidl

–0.992

–1.046

–1.100

–1.261

–1.526

–1.785

–1.000

M22

–0.000

0.002

0.004

0.011

0.023

0.034

0.000

N22

–0.099

–0.105

–0.110

–0.126

–0.153

–0.178

–0.100

σxxin

–0.041

–124.2

–247.9

–615.9

–1219.5

–1811.3

σxxout

–0.183

–58.45

–116.5

–289.2

–572.4

–850.1

Тензорно-векторный способ интерполяции

 y = 0.0; θ = 0.0 σθθin

971.4

971.4

971.4

971.4

971.3

971.3

 σθθout

–977.7

–977.7

–977.7

–977.7

–977.7

–977.7

 M22

–1.624

–1.624

–1.624

–1.624

–1.624

–1.624

 σxxin

291.4

291.4

291.4

291.4

291.4

291.4

 σxxout

–293.3

–293.3

–293.3

–293.3

–293.3

–293.3

 y = 0.5; θ=π2 σθθin

–0.728

–0.728

–0.728

–0.728

–0.728

–0.728

 σθθout

–1.250

–1.250

–1.250

–1.250

–1.250

–1.250

 σθθmidl

–1.002

–1.002

–1.002

–1.002

–1.002

–1.002

–1.000

 M22

–0.001

–0.001

–0.001

–0.001

–0.001

–0.001

0.000

 N22

–0.100

–0.100

–0.100

–0.100

–0.100

–0.100

–0.100

 σxxin

0.214

0.214

0.214

0.214

0.214

0.214

 σxxout

–0.039

–0.039

–0.039

–0.039

–0.039

–0.039

 

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

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

 

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

 

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

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

 

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

Координаты точек, x, м; θ, рад.

Напряжения σ, МПа

Варианты расчета

ПервыйВторой
Сетка узлов
41 × 4146 × 4651 × 5156 × 5641 × 4146 × 4651 × 51

т. A, x = 1.00; θ = 0.00

σθθin

12.36

12.55

12.63

12.64

12.67

12.67

12.66

σθθout

–15.55

–15.53

–15.75

–15.73

–15.74

–15.75

–15.76

т. B, x = 1.00; θ = π

σθθin

49.63

51.44

56.33

57.19

53.56

55.02

56.33

σθθout

–65.20

–67.53

–73.47

–74.32

–71.06

–72.65

–74.08

σxxin

–5.101

–5.183

–5.552

–5.528

–5.574

–5.579

–5.584

σxxout

5.109

5.170

5.549

5.523

5.572

5.578

5.584

т. C, x = 0.00; θ = 0.00

σθθin

13.99

14.05

14.22

14.22

14.23

14.23

14.24

σθθout

–14.15

–14.22

–14.38

–14.38

–14.40

–14.40

–14.40

σxxin

4.286

4.299

4.350

4.349

4.357

4.357

4.358

σxxout

–4.110

–4.100

–4.177

–4.170

–4.188

–4.188

–4.188

т. D, x = 0.00; θ = π

σxxin

3.840

3.891

4.098

4.079

4.111

4.111

4.111

σxxout

–0.690

–0.705

–0.663

–0.664

–0.619

–0.619

–0.618

 

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

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

 

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

 

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

 

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

Координаты точек, x, м; θ, рад.

Напряжения σ, МПа

Величина жесткого смещения в вертикальном направлении, м

0.00

0.10

0.25

0.50

1.0

Первый вариант

т. A, x = 1.00; θ = 0.00

σθθin

12.63

10.31

8.312

6.568

5.035

σθθout

–15.75

–14.41

–13.24

–12.23

–11.34

т. B, x = 1.00; θ = π

σθθin

56.33

45.80

36.72

28.79

21.83

σθθout

–73.47

–58.06

–44.76

–33.16

–22.96

σxxin

–5.552

–4.251

–3.128

–2.149

–1.288

σxxout

5.549

4.328

3.273

2.353

1.544

т. C, x = 0.00; θ = 0.00

σθθin

14.22

12.34

10.71

9.287

8.039

σθθout

–14.38

–12.50

–10.87

–9.447

–8.198

σxxin

4.35

3.953

3.609

3.310

3.046

σxxout

–4.177

–3.253

–2.454

–1.758

–1.146

т. D, x = 0.00; θ = π

σxxin

4.098

3.571

3.117

2.721

2.372

σxxout

–0.663

–1.432

–2.095

–2.6740

–3.183

Второй вариант

т. A, x = 1.00; θ = 0.00

 σθθin

12.66

12.66

12.66

12.66

12.66

 σθθout

–15.76

–15.76

–15.76

–15.76

–15.76

т. B, x = 1.00; θ = π

 σθθin

56.33

56.33

56.33

56.33

56.33

 σθθout

–74.08

–74.08

–74.08

–74.08

–74.08

 σxxin

–5.584

–5.584

–5.584

–5.584

–5.584

 σxxout

5.584

5.584

5.584

5.584

5.584

т. C, x = 0.00; θ = 0.00

 σθθin

4.358

4.358

4.358

4.358

4.358

 σθθout

–4.188

–4.188

–4.188

–4.188

–4.188

 σxxin

14.24

14.24

14.24

14.24

14.24

 σxxout

–14.40

–14.40

–14.40

–14.40

–14.40

т. D, x = 0.00; θ = π

 σxxin

4.111

4.111

4.111

4.111

4.111

 σxxout

–0.618

–0.618

–0.618

–0.618

–0.618

 

Как следует из анализа данных табл. 4, в первом варианте расчета наблюдается существенное изменение параметров НДС, причем фиксируемые изменения пропорциональны величине жесткого смещения оболочки в вертикальном направлении. Так, например, σxxin в точке B уменьшились в 4.3 раза, а σθθout – в 3.2 раза.

Для наглядности полученных вычислений на рис. 5, 6 представлены повариантные графические зависимости изменения значений нормальных напряжений σθθ и σxx в точке B в зависимости от величины жесткого вертикального смещения.

 

Рис. 5. Графики значений нормальных напряжений σθθ в точке B

 

Рис. 6. Графики значений нормальных напряжений σxx в точке B

 

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

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

Заключение

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

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

×

About the authors

M. Yu. Klochkov

Volgograd State Technical University

Author for correspondence.
Email: m.klo4koff@yandex.ru
Russian Federation, Volgograd

V. A. Pshenichkina

Volgograd State Technical University

Email: vap_hm@list.ru
Russian Federation, Volgograd

A. P. Nikolaev

Volgograd State Agricultural University

Email: anpetr40@yandex.ru
Russian Federation, Volgograd

Yu. V. Klochkov

Volgograd State Agricultural University

Email: klotchkov@bk.ru
Russian Federation, Volgograd

O. V. Vakhnina

Volgograd State Agricultural University

Email: ovahnina@bk.ru
Russian Federation, Volgograd

A. S. Andreev

Volgograd State Agricultural University

Email: aandreev.07.1988@gmail.com
Russian Federation, Volgograd

References

  1. Novozhilov V.V. Theory of Thin Shells. St. Petersburg: St. Petersburg Univ. Pub. House, 2010. 378 p. (in Russian)
  2. Balabukh L.I., Alfutov N.A., Usyukin V.I. Structural Mechanics of Rockets. Moscow: Higher School, 1984. 391 p. (in Russian)
  3. Balabukh L.I., Kolesnikov K.S., Zarubin V.S. et al. Fundamentals of Structural Mechanics of Rockets. Moscow: Higher School, 1969. 494 p. (in Russian)
  4. Obraztsov I.F., Vasiliev V.V., Bulychev L.I. et al. Structural mechanics of aircraft. Moscow: Mechanical Engineering, 1986. 536 p. (in Russian)
  5. Postnov V.A., Kharkhurim I.Ya. Finite Element Method in Calculations of Ship Structures. Leningrad: Shipbuilding, 1974. 342 p. (in Russian)
  6. Rickards R.B. Finite Element Method in the Theory of Shells and Plates. Riga: Znatne, 1988. 284 p. (in Russian)
  7. Sekulovich M. Finite Element Method. Moscow: Stroyizdat, 1993. 664 p. (in Russian)
  8. Schöllhammer D., Fries T.P. A Higher-order trace finite element method for shells // Numerical Methods in Engineering, 2021, vol. 122, no. 5. P. 1217–1238.
  9. Yeongbin Ko, Phill-Seung Lee, Klaus-Jürgen Bathe. A new 4-node MITC element for analysis of two-dimensional solids and its formulation in a shell element // Computers & Structures, 2017, vol. 192, pp. 34–49.
  10. Yakupov S.N., Kiyamov H.G., Yakupov N.M. Modeling a synthesized element of complex geometry based upon three-dimensional and two-dimensional finite elements // Lobachevskii J. of Mathematics, 2021, vol. 42, no 9, pp. 2263–2271.
  11. Nguyen Nhung, Waas Anthonym. Nonlinear, finite deformation, finite element analysis // ZAMP. Z. Angew. Math. and Phys, 2016, vol. 67, no. 9, pp. 35/1–35/24.
  12. Gao L., Wang C., Liu Z. et al. Theoretical aspects of selecting repeated unit cell model in micromechanical analysis using displacement-based finite element method // Chinese Journal of Aeronautics, 2017, vol. 30, no. 4, pp. 1417–1426.
  13. Jin He, Jiaxi Zhao, Chenbo Yin. Constitutive equations and stiffness related properties for elastic and hyperelastic solid surfaces: Theories and finite element implementations // International J. of Solids & Structures, 2020, vol. 202, no. 1, pp. 660–671.
  14. Dzhabrailov A.Sh., Nikolaev A.P., Klochkov Yu.V. et al. Finite element algorithm for calculating an ellipsoidal shell taking into account the displacement as a rigid whole // Appl. Math. & Mech., 2022, vol. 86, no. 2, pp. 251–262.
  15. Bakulin V.N. An effective model for layer-by-layer analysis of three-layer irregular shells of rotation of a cylindrical shape // Reports of the Academy of Sci., 2018, vol. 478, no. 2, pp. 148–152.
  16. Bakulin V.N. Model for layer-by-layer analysis of the stress-strain state of three-layer irregular shells of rotation of double curvature // Izv. RAS. MTT, 2020, no. 2, pp. 112–122.
  17. Bakulin V.N. An effective model of load-bearing layers for layer-by-layer analysis of the stress-strain state of three-layer cylindrical irregular shells of revolution // Izv. RAS. MTT, 2020, no. 3, pp. 82–92.
  18. Bakulin V.N. Layer-by-layer analysis of the stress-strain state of irregular three-layer shells of rotation of non-zero Gaussian curvature // Appl. Math. & Mech., 2021, vol. 85, no. 1, pp. 89–105.
  19. Bakulin V.N. Block-layer approach for analyzing the stress-strain state of three-layer irregular cylindrical shells of revolution // Appl. Math. & Mech., 2021, vol. 85, no. 3, pp. 383–395.
  20. Bakulin V.N. Model for analyzing the stress-strain state of three-layer cylindrical shells with rectangular cutouts // Izv. RAS. MTT, 2022, no. 1, pp. 122–132.
  21. Bakulin V.N. Refined model of layer-by-layer analysis of three-layer irregular conical shells // Reports of the Academy of Sciences, 2017, vol. 472, no. 3, pp. 272–277.
  22. Bakulin V.N. Testing a finite element model designed to study the stress-strain state of layered irregular shells // Math. Modeling, 2009, vol. 21, no. 8, pp. 121–128.
  23. Lalin V.V., Rybakov V.A., Ivanov S.S. et al. Mixed finite-element method in V.I. Slivker’s semi-shear thin-walled bar theory // Magazine of Civil Engineering, 2019, vol. 5, no. 89, pp. 79–93.
  24. Klochkov Yu., Pshenichkina V., Nikolaev A. et al. Stress-strain state of elastic shell based on mixed finite element // Magazine of Civil Eng., 2023, vol. 4, no. 120, pp. 12003.
  25. Klochkov Yu.V., Pshenichkina V.A., Nikolaev A.P. et al. Quadrangular finite element in a mixed FEM formulation for the calculation of thin shells of rotation // Structural Mechanics of Engineering Structures and Structures, 2023, vol. 19, no. 1, pp. 64–72.
  26. Magisano D., Liang K., Garcea G. et al. An efficient mixed variational reduced-order model formulation for nonlinear analyses of elastic shells // International J. for Numerical Methods in Engineering, 2018, vol. 113, no. 4, pp. 634–655.
  27. Antonietti P.F., Beirao da Veiga L., Scacchi S. et al. A C1 Virtual element method for the Cahn–Hilliard equation with polygonal meshes // SIAM J. Numer. Anal, 2016, vol. 54, no. 1, pp. 34–56.
  28. Chi H., Talischi C., Lopez-Pamies O. et al. A paradigm for higher order polygonal elements in finite elasticity // Comput. Methods Appl. Mech. Engrg, 2016, vol. 306, pp. 216–251.
  29. Golovanov A.I., Tyuleneva O.N., Shigabutdinov A.F. Finite Element Method in the Statics and Dynamics of Thin-Walled Structures. Moscow: Fizmatlit, 2006. 392 p. (in Russian)
  30. Skopinsky V.N. Stresses in Intersecting Shells. Moscow: Fizmatlit, 2008. 400 p. (in Russian)
  31. Klochkov Yu.V., Nikolaev A.P., Sobolevskaya T.A. et al. The calculation of the ellipsoidal shell-based FEM with vector interpolation of displacements when the variable parameterisation of the middle surface // Lobachevskii Journal of Mathematics, 2020, vol. 41, no. 3, pp. 373–381.
  32. Dzhabrailov A.Sh., Nikolaev A.P., Klochkov Yu.V. et al. Taking into account displacement as a rigid body in the FEM algorithm when calculating shells of revolution // Izv. RAS. MTT, 2023, no. 6, pp. 23–38.

Supplementary files

Supplementary Files
Action
1. JATS XML
2. Fig. 1. Calculation diagram of an elliptical ring with hinged support

Download (22KB)
3. Fig. 2. Calculation diagram of elliptical ring with spring support

Download (23KB)
4. Fig. 3. Calculation diagram of an elliptical cylinder with hinged support

Download (29KB)
5. Fig. 4. Calculation diagram of elliptical cylinder with spring support

Download (29KB)
6. Fig. 5. Plots of normal stress values σθθ at point B

Download (26KB)
7. Fig. 6. Plots of normal stress values σxx at point B

Download (27KB)

Copyright (c) 2024 Russian Academy of Sciences

Согласие на обработку персональных данных с помощью сервиса «Яндекс.Метрика»

1. Я (далее – «Пользователь» или «Субъект персональных данных»), осуществляя использование сайта https://journals.rcsi.science/ (далее – «Сайт»), подтверждая свою полную дееспособность даю согласие на обработку персональных данных с использованием средств автоматизации Оператору - федеральному государственному бюджетному учреждению «Российский центр научной информации» (РЦНИ), далее – «Оператор», расположенному по адресу: 119991, г. Москва, Ленинский просп., д.32А, со следующими условиями.

2. Категории обрабатываемых данных: файлы «cookies» (куки-файлы). Файлы «cookie» – это небольшой текстовый файл, который веб-сервер может хранить в браузере Пользователя. Данные файлы веб-сервер загружает на устройство Пользователя при посещении им Сайта. При каждом следующем посещении Пользователем Сайта «cookie» файлы отправляются на Сайт Оператора. Данные файлы позволяют Сайту распознавать устройство Пользователя. Содержимое такого файла может как относиться, так и не относиться к персональным данным, в зависимости от того, содержит ли такой файл персональные данные или содержит обезличенные технические данные.

3. Цель обработки персональных данных: анализ пользовательской активности с помощью сервиса «Яндекс.Метрика».

4. Категории субъектов персональных данных: все Пользователи Сайта, которые дали согласие на обработку файлов «cookie».

5. Способы обработки: сбор, запись, систематизация, накопление, хранение, уточнение (обновление, изменение), извлечение, использование, передача (доступ, предоставление), блокирование, удаление, уничтожение персональных данных.

6. Срок обработки и хранения: до получения от Субъекта персональных данных требования о прекращении обработки/отзыва согласия.

7. Способ отзыва: заявление об отзыве в письменном виде путём его направления на адрес электронной почты Оператора: info@rcsi.science или путем письменного обращения по юридическому адресу: 119991, г. Москва, Ленинский просп., д.32А

8. Субъект персональных данных вправе запретить своему оборудованию прием этих данных или ограничить прием этих данных. При отказе от получения таких данных или при ограничении приема данных некоторые функции Сайта могут работать некорректно. Субъект персональных данных обязуется сам настроить свое оборудование таким способом, чтобы оно обеспечивало адекватный его желаниям режим работы и уровень защиты данных файлов «cookie», Оператор не предоставляет технологических и правовых консультаций на темы подобного характера.

9. Порядок уничтожения персональных данных при достижении цели их обработки или при наступлении иных законных оснований определяется Оператором в соответствии с законодательством Российской Федерации.

10. Я согласен/согласна квалифицировать в качестве своей простой электронной подписи под настоящим Согласием и под Политикой обработки персональных данных выполнение мною следующего действия на сайте: https://journals.rcsi.science/ нажатие мною на интерфейсе с текстом: «Сайт использует сервис «Яндекс.Метрика» (который использует файлы «cookie») на элемент с текстом «Принять и продолжить».