Three-Field FEM in Shell Calculations with Options for Interpolation of the Sought Values
- 作者: Klochkov M.Y.1, Pshenichkina V.A.1, Nikolaev A.P.2, Klochkov Y.V.2, Vakhnina O.V.2, Andreev A.S.2
-
隶属关系:
- Volgograd State Technical University
- Volgograd State Agricultural University
- 期: 卷 88, 编号 5 (2024)
- 页面: 797-820
- 栏目: Articles
- URL: https://journals.rcsi.science/0032-8235/article/view/280969
- DOI: https://doi.org/10.31857/S0032823524050109
- EDN: https://elibrary.ru/JPFBTH
- ID: 280969
如何引用文章
全文:
详细
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.
全文:
Введение
Оболочки различных конфигураций в последнее время получают все более широкое распространение во многих инженерных сооружениях (трубопроводы, резервуары, покрытия, перекрытия, летательные аппараты, архитектурные формы и другие).
Теория расчета тонких оболочек в настоящее время является достаточно развитой [1–4] при различных видах нагружения. Аналитические решения уравнений теории тонких оболочек оказались возможными лишь в некоторых частных случаях, далеких от инженерной практики. Поэтому актуальными оказались разработки численных методов решения уравнений теории оболочек. Среди численных методов решения задач деформирования инженерных конструкций получил развитие и широкое использование метод конечных элементов (МКЭ) [5–7]. Наиболее широкое использование при расчете оболочек получил МКЭ в формулировке метода перемещений [8–14]. МКЭ в формулировке метода перемещений использовался в расчетах трехслойных оболочек [15–22]. Для обоснования разработанной методики расчета трехслойных оболочек в [21, 22] выполнено сравнение полученных результатов с аналитическими решениями других исследователей.
В расчетах оболочек МКЭ широко использовался и в смешанной формулировке [23–26]. В последние годы разрабатывается виртуальный МКЭ с объемными конечными элементами на пространственных сетках дискретизации [27, 28].
При использовании МКЭ в формулировке метода перемещений применялись в конечных элементах кинематические узловые неизвестные. При смешанной формулировке МКЭ в конечных элементах узловыми неизвестными принимались перемещения и напряжения с выполнением условий совместности между конечными элементами в сетке дискретизации.
При проектировании тонкостенных оболочек необходимо знать численные значения компонент векторов перемещений, а также компонент тензоров напряжений. Достижению этой цели служит развитие численных методов расчета, в частности МКЭ в виде трехпольной формулировки, которая позволяет без дополнительных вычислительных процедур одновременно получать численные значения компонент вектора перемещения, изгибающих моментов, продольных сил, деформаций и искривлений срединной поверхности рассчитываемой оболочки. При этом следует применять вычислительные технологии, позволяющие учитывать смещение конечного элемента как твердого целого.
Получение аппроксимирующих выражений искомых величин с учетом смещения как твердого тела для конечного элемента в форме цилиндрической панели показано в [29, 30], для трехслойных оболочек в [15, 16].
Проблеме учета жестких смещений в МКЭ в векторной формулировке посвящены работы [31, 32].
С целью решения проблем численного анализа НДС оболочек в настоящей работе излагается разработанный трехпольный вариант МКЭ со следующими полями узловых неизвестных:
- перемещения и их первые производные;
- деформации и искривления срединной поверхности;
- усилия и моменты срединной поверхности.
Аппроксимация искомых величин использована в двух вариантах.
В первом варианте традиционные аппроксимирующие выражения использовались для аппроксимации непосредственно компонент векторов перемещений и компонент тензоров деформаций, искривлений, усилий и моментов.
Во втором варианте традиционные аппроксимирующие процедуры применялись к искомым векторам и тензорам, и после координатных преобразований в используемой криволинейной системе координат определялись аппроксимирующие выражения для компонент искомых векторов и тензоров.
Используемая векторно-тензорная форма интерполяционной процедуры позволяет эффективно решить проблему учета смещений оболочек как абсолютно твердых тел.
Геометрия оболочки
Параметризация срединной поверхности оболочки является важным аспектом процесса разработки вычислительного алгоритма расчета тонкостенных объектов. В представленном исследовании в качестве криволинейных координат срединной поверхности были использованы осевая координата x и полярный угол θ. Таким образом, радиус-вектор срединной поверхности оболочки может быть задан векторной функцией следующего вида
, (2.1)
где θ – полярный угол, отсчитываемый в плоскости yOz от оси Oz против хода часовой стрелки; r = r (x, θ) – функция, зависящая от конкретного типа поверхности.
Так, например, для оболочки в форме эллипсоида данная функция имеет вид
, (2.2)
где a, b, c – параметры эллипсоида.
Для эллиптического цилиндра эта формула зависит только от полярного угла
(2.3)
Базисные векторы точки M0 срединной поверхности оболочки в исходном состоянии можно получить стандартным образом с использованием формул дифференциальной геометрии
, (2.4)
где .
Векторы базиса точки M0ζ, расположенной в произвольном слое оболочки на расстоянии ζ от точки M0 в исходном состоянии определяются по формулам
, (2.5)
где .
В процессе деформирования оболочки под действием приложенной внешней нагрузки точки M0 и M0ζ займут новые положения M и Mζ, определяемые векторами перемещений и соответственно
(2.6)
Входящий в структуру формулы вектора орт нормали к деформированной срединной поверхности определяется векторным произведением
, (2.7)
где .
Векторы базиса в точке Mζ деформированного состояния оболочки определяются зависимостями
(2.8)
Метрические тензоры точек M0ζ и Mζ определяются компонентами
, (2.9)
где нижние индексы α и β последовательно принимают значения 1, 2.
Ковариантные компоненты тензора деформаций в точке Mζ могут быть получены посредством использования соотношения механики сплошной среды
. (2.10)
При использовании гипотезы тонких оболочек Кирхгофа–Лява [1] соотношения (2.10) могут быть представлены суммой
, (2.11)
где – деформации и искривления в точках Mζ и M оболочки.
Дискретная модель объекта
Для построения дискретной модели оболочки был выбран четырехузловой конечный элемент в виде фрагмента срединной поверхности с узлами i, j, k, l.
Криволинейные координаты x и θ дискретного элемента определяются с использованием билинейных функций через координаты узлов посредством матричных соотношений
, (3.1)
где .
В качестве искомых кинематических и силовых параметров конечного элемента выбираются следующие столбцы узловых неизвестных: – строки узловых перемещений, каждая из которых имеет вид , где под q понимается компонента вектора перемещения v1, v2 или v; – строки узловых деформаций и искривлений срединной поверхности , где под λαβ понимается деформация εαβ или искривления ℵαβ срединной поверхности; – строки узловых продольных сил и изгибающих моментов . Здесь под Tαβ понимается контравариантная компонента тензора внутренних усилий или тензора моментов.
В разработанных к настоящему времени вычислительных комплексах (ANSYS, NASTRAN, ABAQUS и других) и конечно-элементных алгоритмах традиционно используется покомпонентная интерполяционная процедура, отличительной особенностью которой является применение интерполяционного выражения отдельно к каждой компоненте вектора перемещения или к каждой компоненте тензорных величин, таких как тензор деформаций, искривлений, сил и моментов. Согласно этой интерполяционной процедуре, можно записать следующие интерполяционные выражения
, (3.2)
где – матрица-строка, содержащая произведения полиномов Эрмита третьей степени [24, 25].
При использовании билинейных функций формы для аппроксимации деформаций условия неразрывности деформаций не выполняются из-за наличия деформаций сдвига εxy в узле k конечного элемента.
С использованием интерполяционных зависимостей (3.2) формируются аппроксимирующие матричные выражения
, (3.3)
где – строка деформаций и искривлений внутренней точки конечного элемента; – строка усилий и моментов внутренней точки конечного элемента; – строка деформаций и искривлений внутренней точки, определяемых функциями перемещений.
Как следует из (3.2), каждая из компонент вектора и тензора интерполируется через узловые значения этой же самой компоненты. Данная аппроксимация является корректной, если векторы и тензоры определены в декартовой системе координат.
Если векторы и тензоры определены в криволинейной системе координат, то необходимо использовать тензорно-векторную форму интерполяционной процедуры, при которой выполняется интерполирование не отдельных компонентов вектора или отдельных компонент тензоров деформаций, искривлений, усилий или моментов, а непосредственно самого вектора перемещения и вышеупомянутых тензоров посредством узловых значений векторов и тензоров.
При реализации трехпольного варианта МКЭ применительно к анализу НДС оболочек тензорно-векторная форма интерполяционной процедуры предусматривает использование следующих интерполяционных выражений
, (3.4)
где – строки узловых тензоров деформаций и искривлений срединной поверхности; ; – строки узловых тензоров продольных сил и изгибающих моментов; – строка узловых векторов перемещения и их производных первого порядка.
Тензоры деформаций и искривлений, усилий и моментов, входящих в (3.4), можно представить компонентами диадных произведений соответствующих векторов локальных базисов, например
(3.5)
где
матрица имеет вид , но ее компоненты, отличные от нуля, являются диадными произведениями контравариантных векторов базисов узловых точек;
Векторы базисов внутренних и узловых точек конечного элемента выражаются через орты декартовой системы координат посредством следующих матричных соотношений
, (3.6)
где , , . Верхний индекс ρ обозначает конкретный из узлов КЭ i, j, k, l.
В результате обращения первого соотношения (3.6) и подстановки его во второе соотношение (3.6) базисные векторы узлов конечного элемента выражаются через базисные векторы его внутренней точки
(3.7)
С учетом (3.7) диадные произведения узлов КЭ также могут быть выражены через диадные произведения точки, принадлежащей внутренней области конечного элемента
, (3.8)
где .
Принимая во внимание (3.8), матрица примет вид
(3.9)
При учете (3.5)–(3.9) соотношения (3.4) могут быть представлены выражениями
(3.10)
Из (3.10) можно получить необходимые интерполяционные зависимости для компонент тензоров деформаций, искривлений, продольных сил и изгибающих моментов, например
(3.11)
Получение интерполяционных соотношений для компонент вектора перемещения в соответствие с векторной формой интерполяционной процедуры изложено [24, 25].
На основе интерполяционных зависимостей (3.11) формируются аппроксимирующие выражения
. (3.12)
Анализируя (3.11), можно констатировать, что отдельная компонента любого тензора (деформаций, искривлений, усилий и моментов) в точке, принадлежащей внутренней области используемого элемента дискретизации, зависит от узловых значений всех его компонент, то есть от полного набора узловых варьируемых параметров соответствующего тензора.
При общепринятой форме интерполяционной процедуры (3.2) отдельная компонента тензора, например, компонента тензора искривлений ℵ11 является функцией узловых значений только этой же самой компоненты, т.е. и не зависит от узловых значений остальных компонент данного тензора .
Также особенностью новых интерполяционных зависимостей (3.11) является присутствие в них параметров применяемой в конкретном расчете оболочки криволинейной системы координат, что не наблюдается в традиционных интерполяционных соотношениях (3.2).
Матрица жесткости конечного элемента
Для получения матрицы жесткости конечного элемента используется функционал [24, 25], основанный на равенстве работ заданных сил на перемещениях и внутренних усилий на деформациях и искривлениях, дополненный условием равенства нулю работы на деформациях и искривлениях невязки, определяемой как разность между внутренними усилиями, вводимыми по определению, и внутренними усилиями, выражаемыми через деформации и искривления
, (4.1)
где , – матрицы-строки продольных сил и изгибающих моментов, а также деформационных параметров в точке срединной поверхности оболочки; {εℵk} – столбец, элементы которого определяются соотношениями Коши; ; – матрицы-строки компонент вектора перемещения и вектора внешней поверхностной нагрузки.
Входящая в (4.1) матрица [C] определяет связь между столбцами {NM} и {εℵ}
(4.2)
При учете аппроксимирующих выражений для кинематических {U}, {εℵ} и силовых {NM} искомых неизвестных в зависимости от используемого варианта аппроксимации (3.2) или (3.11) функционал (4.1) с учетом (4.2) может быть преобразован к виду
(4.3)
где матрицы , и компонуются на основе интерполяционных выражений (3.3) или (3.12), нижний индекс γ принимает значения 1, 2, причем γ = 1 соответствует использованию общепринятой покомпонентной интерполяционной процедуре (3.2), а γ = 2 соответствует применению разработанной тензорно-векторной форме интерполяции искомых величин (3.11) в трехпольном варианте МКЭ.
Входящая в (4.3) матрица определяет связь между столбцами кинематических узловых неизвестных в локальной и глобальной системах координат.
Осуществляя минимизацию (4.3) по кинематическим , {εℵy}T и силовым {NMy}T искомым неизвестным, можно записать систему матричных уравнений
(4.4)
или в более компактном виде
(4.5)
где .
Выражая из третьего и второго уравнений (4.5) столбцы {εℵy}T и {NMy}T, можно получить следующие матричные зависимости
(4.6)
Осуществляя подстановку (4.6) в первое уравнение (4.5), можно записать следующее матричное выражение
, (4.7)
где – искомая матрица жесткости размером 36 × 36 конечного элемента, используемого для построения дискретной модели оболочки.
Глобальная матрица жесткости всей исследуемой оболочки, представляющая собой матрицу системы линейных алгебраических уравнений (СЛАУ), компонуется из отдельных матриц жесткостей [K] посредством использования матрицы индексов [5].
Полученные в результате решения методом Гаусса глобальной СЛАУ столбцы посредством (4.6), позволяют вычислять требуемые для анализа НДС оболочки столбцы продольных сил и изгибающих моментов {εℵy}, а также столбцы деформаций и искривлений {NMy} в любой интересующей исследователя узловой точке рассчитываемой конструкции.
Вычислительные эксперименты
Эксперимент 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 даются значения и , определенные соотношениями
Таблица 1. Значения нормальных напряжений, изгибающих моментов и продольных сил в точке приложения нагрузки и в точке шарнирного опирания
Координаты точек, y, м; θ, рад. | σ, 10–2 МПа; M, Н · см; N, Н | Способ интерполяции | Значения M22, N22, из условия равновесия | |||||||
Общепринятый | Тензорно-векторный | |||||||||
Сетка узлов | ||||||||||
121 × 2 | 151 × 2 | 181 × 2 | 211 × 2 | 121 × 2 | 151 × 2 | 181 × 2 | 211 × 2 | |||
y = 0.0; θ = 0.0 | 971.4 | 971.4 | 971.4 | 971.4 | 971.3 | 971.3 | 971.4 | 971.4 | – | |
–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 | – | |
291.4 | 291.4 | 291.4 | 291.4 | 291.4 | 291.4 | 291.4 | 291.4 | – | ||
–293.3 | –293.3 | –293.3 | –293.3 | –293.3 | –293.3 | –293.3 | –293.3 | – | ||
y = 0.5; | 0.319 | –0.405 | –0.703 | –0.840 | 0.949 | –0.091 | –0.524 | –0.728 | – | |
–2.151 | –1.508 | –1.247 | –1.129 | –2.834 | –1.849 | –1.441 | –1.250 | – | ||
–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 | |
1.825 | 0.612 | 0.157 | –0.041 | 2.135 | 0.993 | 0.475 | 0.214 | – | ||
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, из условия равновесия | |||||
0.00 | 0.05 | 0.10 | 0.25 | 0.50 | 0.75 | |||
Общепринятый способ интерполяции | ||||||||
y = 0.0; θ = 0.0 | 971.4 | 971.2 | 971.1 | 970.6 | 969.8 | 969.0 | – | |
–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 | – | |
291.4 | 291.4 | 291.35 | 291.25 | 291.1 | 290.9 | – | ||
–293.3 | –293.3 | –293.2 | –293.0 | –292.7 | –292.4 | – | ||
y = 0.5; | –0.840 | –0.246 | 0.346 | 2.107 | 4.994 | 7.825 | – | |
–1.129 | –0.770 | –0.414 | 0.648 | 2.390 | 4.098 | – | ||
–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 | |
–0.041 | –124.2 | –247.9 | –615.9 | –1219.5 | –1811.3 | – | ||
–0.183 | –58.45 | –116.5 | –289.2 | –572.4 | –850.1 | – | ||
Тензорно-векторный способ интерполяции | ||||||||
y = 0.0; θ = 0.0 | 971.4 | 971.4 | 971.4 | 971.4 | 971.3 | 971.3 | – | |
–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 | – | |
291.4 | 291.4 | 291.4 | 291.4 | 291.4 | 291.4 | – | ||
–293.3 | –293.3 | –293.3 | –293.3 | –293.3 | –293.3 | – | ||
y = 0.5; | –0.728 | –0.728 | –0.728 | –0.728 | –0.728 | –0.728 | – | |
–1.250 | –1.250 | –1.250 | –1.250 | –1.250 | –1.250 | – | ||
–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 | |
0.214 | 0.214 | 0.214 | 0.214 | 0.214 | 0.214 | – | ||
–0.039 | –0.039 | –0.039 | –0.039 | –0.039 | –0.039 | – |
Анализ данных табл. 2 показывает, что в первом варианте расчета параметры НДС оболочки в точке пружинного опирания претерпевают весьма значительные изменения, которые возрастают по мере увеличения смещения оболочки как абсолютно твердого тела. Так, при смещении оболочки на величину 0.10 м напряжения уже меняют свой знак, а σ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 × 41 | 46 × 46 | 51 × 51 | 56 × 56 | 41 × 41 | 46 × 46 | 51 × 51 | |||
т. A, x = 1.00; θ = 0.00 | 12.36 | 12.55 | 12.63 | 12.64 | 12.67 | 12.67 | 12.66 | ||
–15.55 | –15.53 | –15.75 | –15.73 | –15.74 | –15.75 | –15.76 | |||
т. B, x = 1.00; θ = π | 49.63 | 51.44 | 56.33 | 57.19 | 53.56 | 55.02 | 56.33 | ||
–65.20 | –67.53 | –73.47 | –74.32 | –71.06 | –72.65 | –74.08 | |||
–5.101 | –5.183 | –5.552 | –5.528 | –5.574 | –5.579 | –5.584 | |||
5.109 | 5.170 | 5.549 | 5.523 | 5.572 | 5.578 | 5.584 | |||
т. C, x = 0.00; θ = 0.00 | 13.99 | 14.05 | 14.22 | 14.22 | 14.23 | 14.23 | 14.24 | ||
–14.15 | –14.22 | –14.38 | –14.38 | –14.40 | –14.40 | –14.40 | |||
4.286 | 4.299 | 4.350 | 4.349 | 4.357 | 4.357 | 4.358 | |||
–4.110 | –4.100 | –4.177 | –4.170 | –4.188 | –4.188 | –4.188 | |||
т. D, x = 0.00; θ = π | 3.840 | 3.891 | 4.098 | 4.079 | 4.111 | 4.111 | 4.111 | ||
–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 | 12.63 | 10.31 | 8.312 | 6.568 | 5.035 | |
–15.75 | –14.41 | –13.24 | –12.23 | –11.34 | ||
т. B, x = 1.00; θ = π | 56.33 | 45.80 | 36.72 | 28.79 | 21.83 | |
–73.47 | –58.06 | –44.76 | –33.16 | –22.96 | ||
–5.552 | –4.251 | –3.128 | –2.149 | –1.288 | ||
5.549 | 4.328 | 3.273 | 2.353 | 1.544 | ||
т. C, x = 0.00; θ = 0.00 | 14.22 | 12.34 | 10.71 | 9.287 | 8.039 | |
–14.38 | –12.50 | –10.87 | –9.447 | –8.198 | ||
4.35 | 3.953 | 3.609 | 3.310 | 3.046 | ||
–4.177 | –3.253 | –2.454 | –1.758 | –1.146 | ||
т. D, x = 0.00; θ = π | 4.098 | 3.571 | 3.117 | 2.721 | 2.372 | |
–0.663 | –1.432 | –2.095 | –2.6740 | –3.183 | ||
Второй вариант | ||||||
т. A, x = 1.00; θ = 0.00 | 12.66 | 12.66 | 12.66 | 12.66 | 12.66 | |
–15.76 | –15.76 | –15.76 | –15.76 | –15.76 | ||
т. B, x = 1.00; θ = π | 56.33 | 56.33 | 56.33 | 56.33 | 56.33 | |
–74.08 | –74.08 | –74.08 | –74.08 | –74.08 | ||
–5.584 | –5.584 | –5.584 | –5.584 | –5.584 | ||
5.584 | 5.584 | 5.584 | 5.584 | 5.584 | ||
т. C, x = 0.00; θ = 0.00 | 4.358 | 4.358 | 4.358 | 4.358 | 4.358 | |
–4.188 | –4.188 | –4.188 | –4.188 | –4.188 | ||
14.24 | 14.24 | 14.24 | 14.24 | 14.24 | ||
–14.40 | –14.40 | –14.40 | –14.40 | –14.40 | ||
т. D, x = 0.00; θ = π | 4.111 | 4.111 | 4.111 | 4.111 | 4.111 | |
–0.618 | –0.618 | –0.618 | –0.618 | –0.618 |
Как следует из анализа данных табл. 4, в первом варианте расчета наблюдается существенное изменение параметров НДС, причем фиксируемые изменения пропорциональны величине жесткого смещения оболочки в вертикальном направлении. Так, например, в точке B уменьшились в 4.3 раза, а – в 3.2 раза.
Для наглядности полученных вычислений на рис. 5, 6 представлены повариантные графические зависимости изменения значений нормальных напряжений σθθ и σxx в точке B в зависимости от величины жесткого вертикального смещения.
Рис. 5. Графики значений нормальных напряжений σθθ в точке B
Рис. 6. Графики значений нормальных напряжений σxx в точке B
На графиках отчетливо прослеживается тенденция к увеличению погрешности вычисления напряжений σθθ и σxx в первом варианте расчета по мере увеличения жесткого вертикального смещения оболочки. Во втором варианте расчета параметры НДС остаются неизменными, чему свидетельствуют строго горизонтальные графики σθθ и σxx.
Из анализа данных табл. 4 видно, что во втором варианте расчета параметры НДС оставались абсолютно неизменными при любом значении величины вертикального жесткого смещения, что приводит к выводу о том, что разработанная тензорно-векторная форма интерполяционной процедуры искомых величин в трехпольном варианте МКЭ позволяет полностью учитывать смещения исследуемых оболочек как абсолютно твердых тел.
Заключение
Разработан векторно-тензорный способ интерполяционной процедуры для трехпольного варианта МКЭ, позволяющий выразить компоненты силовых (продольных сил и изгибающих моментов) и кинематических (деформаций и искривлений срединной поверхности) искомых неизвестных через полный набор узловых силовых или кинематических параметров.
Выполненные численные эксперименты по расчету оболочек с эллиптическим поперечным сечением, имеющих возможность под действием заданной нагрузки смещаться как абсолютно твердые тела, показали, что применение разработанной тензорно-векторной формы интерполяционной процедуры позволяет корректным образом учитывать жесткие смещения оболочек и вычислять искомые параметры НДС без какой-либо погрешности.
作者简介
M. Klochkov
Volgograd State Technical University
编辑信件的主要联系方式.
Email: m.klo4koff@yandex.ru
俄罗斯联邦, Volgograd
V. Pshenichkina
Volgograd State Technical University
Email: vap_hm@list.ru
俄罗斯联邦, Volgograd
A. Nikolaev
Volgograd State Agricultural University
Email: anpetr40@yandex.ru
俄罗斯联邦, Volgograd
Yu. Klochkov
Volgograd State Agricultural University
Email: klotchkov@bk.ru
俄罗斯联邦, Volgograd
O. Vakhnina
Volgograd State Agricultural University
Email: ovahnina@bk.ru
俄罗斯联邦, Volgograd
A. Andreev
Volgograd State Agricultural University
Email: aandreev.07.1988@gmail.com
俄罗斯联邦, Volgograd
参考
- Novozhilov V.V. Theory of Thin Shells. St. Petersburg: St. Petersburg Univ. Pub. House, 2010. 378 p. (in Russian)
- Balabukh L.I., Alfutov N.A., Usyukin V.I. Structural Mechanics of Rockets. Moscow: Higher School, 1984. 391 p. (in Russian)
- 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)
- Obraztsov I.F., Vasiliev V.V., Bulychev L.I. et al. Structural mechanics of aircraft. Moscow: Mechanical Engineering, 1986. 536 p. (in Russian)
- Postnov V.A., Kharkhurim I.Ya. Finite Element Method in Calculations of Ship Structures. Leningrad: Shipbuilding, 1974. 342 p. (in Russian)
- Rickards R.B. Finite Element Method in the Theory of Shells and Plates. Riga: Znatne, 1988. 284 p. (in Russian)
- Sekulovich M. Finite Element Method. Moscow: Stroyizdat, 1993. 664 p. (in Russian)
- 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.
- 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.
- 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.
- 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.
- 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.
- 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.
- 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.
- 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.
- 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.
- 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.
- 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.
- 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.
- 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.
- 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.
- 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.
- 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.
- 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.
- 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.
- 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.
- 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.
- 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.
- 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)
- Skopinsky V.N. Stresses in Intersecting Shells. Moscow: Fizmatlit, 2008. 400 p. (in Russian)
- 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.
- 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.
