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

Во введении обоснована актуальность темы исследования, указаны цели
изадачиисследования,научнаяновизна,практическаязначимость
диссертационной работы, применяемые методы исследований, методология
исследования, положения, выносимые на защиту.
В первой главе приводится краткий обзор существующих в настоящее
время отечественных и зарубежных работ по исследуемой теме. В большинстве
случаев авторы используют в своих работах классическую теорию пластин и
оболочек, которая была проверена временем и нашла широкое применение на
практике. При использовании классической теории возникает погрешность,
которая связанна с гипотезой Кирхгофа. Эту проблему можно решить, если
будет проведён учет деформаций поперечного сдвига. Проведённый анализ
работ также показал, что при расчете инженерных конструкций, у которых
наблюдается значительный прогиб элементов, решение в линейной постановке
не вполне корректно, так как приводит к погрешности результатов. Решение
такого рода задач в геометрически нелинейной постановке с учетом сдвиговой
теории приводит к более корректному результату. В работах достаточно полно
изучен вопрос с решением задач в геометрически линейной постановке. Расчёт
оболочек с учётом геометрической и физической нелинейности на основе
методаконечныхэлементовтребуетдальнейшегоразвитияи
совершенствования.
Применение стандартной интерполяционной процедурыотдельных
компонент шагового вектора перемещения и шагового вектора угла наклона
нормали при расчетах в криволинейных системах координат приводит к
погрешностям вычислений. В таких случаях предпочтительнее использование
векторной формы интерполяции шаговых искомых величин, так как она
автоматически учитывает смещение конечного элемента как жесткого целого,
что приводит к получению более корректных результатов расчёта.
Во второй главе, посвященной анализу НДС тонких оболочек в
линейной постановке с учётом сдвиговых деформаций, приводится вывод
основных геометрических и физических соотношений упругих оболочек.
Представлен алгоритм формирования матрицы жесткости четырёхузлового
конечногоэлементаприиспользованиистандартнойдляМКЭ
интерполяционнойпроцедурыивекторной
альтернативнойформы

интерполяции искомых величин. Отсчет вектора наклона нормали 
осуществляется от исходного состояния нормали. На тестовых примерах
расчёта показана эффективность применения векторной интерполяционной
процедуры при использовании четырехузлового конечного элемента с
узловыми неизвестными в виде компонент вектора перемещения и их
производных, а также компонент вектора угла наклона нормали.
Пример 1. В качестве примера решения задачи в линейной постановке
было исследовано НДС цилиндра, загруженного вдоль образующей равномерно
распределенной нагрузкой интенсивности q и имеющего по торцам шарнирные
опоры, препятствующие вертикальному смещению (рисунок 1), которые в по-
следствиибылизамененына
пружинные. Торцы цилиндра имели
скользящую заделку. Были приняты
следующие исходные данные: q = 1
H/см; L = 1,5 м; R = 0,4 м; t = 0,6 см;
E = 2·107 Н/см2; v = 0,3.
Рисунок 1 – Расчетная схема цилиндра с
шарнирными опорами (с пружинными опорами)Расчеты были выполнены в двух
вариантах: при реализации стандартной и векторной форм интерполяции
искомых неизвестных. Результаты повариантных расчетов представлены в
таблице 1, из анализа которой следует, что при шарнирном опирании оба
варианта интерполяции искомых величин приводят к примерно одинаковым
результатам в численных значениях расчетных величин.
Таблица 1 – Результаты расчётов с использованием стандартной и векторной
интерполяционных процедур при различных величинах жестких смещений цилиндрической
оболочки
Интерполяционная процедура
Напряжения,
(х, см; t, рад)
Координаты

СтандартнаяВекторная
Н/см2
точек

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

0.0200.0500.00.0200.0500.0
(шарнирное опирание)(шарнирное опирание)
 11в111.5189.7301.8111.0111.0111.0
 22в
Т. 1 (20.0; 0.0)

384.2634.6994.0382.3382.3382.3
 13в-32.2-31.1-32.0-32.3-32.3-32.3
н
-108.8-30.881.1-108.5-108.5-108.5
 22н-382.0-131.6227.7-380.7-380.7-380.7
 13н-32.2-32.1-32.0-32.3-32.3-32.3
в
460.6374.3250.4455.5455.5455.5
 22в717.8454.777.0696.9696.9696.9
Т. 3 (20.0; π)

в
-253.2-248.8-242.5-254.8-254.8-254.8
 11н-457.6-527.4-627.7-461.4-461.4-461.4
 22н-826.3-1062.1-1400.4-842.1-842.1-842.1
н
-253.2-248.8-242.5-254.8-254.8-254.8

В случае замены шарнирных опор на пружинные (рисунок 1) оболочка
получает возможность смещаться в вертикальном направлении как абсолютно
жесткое тело. Изменяя жесткость пружины, можно получить различные
величины жесткого смещения.
Анализируя данные таблицы 1, можно убедиться, что при использовании
стандартнойинтерполяционнойпроцедурыпараметрыНДСцилиндра
существенно меняются в зависимости от величины жесткого смещения.
Причем, с увеличением жесткого смещения значения напряжений возрастают и
достигают неприемлемых значений, некоторые из них меняют свой знак. При
векторной форме интерполяционной процедуры параметры НДС цилиндра не
изменяются, несмотря на существенную величину жесткого смещения.
Полученный результат объясняется тем, что использование векторной формы
интерполяционной процедуры позволяет автоматически учитывать смещения
КЭ как жесткого целого при использовании криволинейных систем координат.
В третьей главе получены основные геометрические соотношения
тонких оболочек с учётом сдвиговых деформаций в геометрически нелинейной
постановке, определяющие связь между приращениями деформаций и компо-
нентамишаговоговектора
перемещения,компонентами
шаговоговектораугланаклона
нормали,атакжепервыми
производнымивышеуказанных
компонент. Основные соотношения
Рисунок 2 – Перемещение точек срединной
были получены на основе анализа
поверхности М0 и точки произвольного слоя
М0ζ при шаговом нагружениитрех состояний системы: исходного
и двух деформированных за j шагов и на ( j  1) шаге нагружения (рисунок 2).
Суммарный и шаговый векторы перемещений, а также их первые
производныепоестественнымкриволинейнымкоординатамs1иs2
определяются компонентами, отнесенными к базису исходного состояния
 
v  ue10  ve20  wen0 ; w  ue10  ve20  wen0 ;
 
v,  1 e10  2 e20  n en0 ; w,  l1 e10  l2 e20  ln en0 .(1)

где   s1 , s2 .

Радиус-векторы R 0 , R  и R * , определяющие соответственно положения
точек M 0 , M  и M * (рисунок 2) могут быть представлены выражениями
 
R 0  R 0  en0 ; R   R 0  V ; R *  R   W . (2)

Входящие в (2) суммарный V и шаговый W векторы перемещений точек
M 0 и M  можно записать в виде сумм
    
V  v    v   ( 1e10   2e20 ); W  w    w   ( 1e10   2 e20 ).(3)
Приращения деформаций в точке произвольного слоя оболочки,
отстоящего от срединной поверхности на расстоянии  , находятся с помощью
уравнений механики сплошной среды
   g 
*
 g   / 2,(4)
  
где g
*
 g*  g * ; g  g  g  .

Cтолбцы узловых варьируемых параметров четырёхузлового конечного
элемента в локальной –  , и глобальной s1 , s 2 системах координат можно
записать в следующем виде

U л T
 u yл  v yл  w yл   1y   y2  ;
TTTTT 

 112
y
1 441121121 41 4

U ГT
 u yГ  v yГ  w yГ   1y   y2  .
TTTTT 
(5)
 112
y
1 441121121 41 4

Искомые величиныU внутренней точки конечного элемента
51

аппроксимируются через узловые неизвестные с помощью матричного
выражения
U   AU ,л
y
(6)
51544441

гдеA- матрица, элементами которой являются произведения
544

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

ФL          k   dV   U   P k PdF .
TT
(7)
V15 5151 F15 5151 

Коэффициент k в (7) принимает значения 1 или 1/2, причём при k  1
соответствует равенству возможных, а при k  1 / 2 – равенству действительных
работ внешних и внутренних сил на ( j  1) -м шаге нагружения.
Пример 2. Была решена задача
поопределениюНДСжестко
защемленнойпообразующим
цилиндрическойпанели,
загруженнойвцентрепролета
сосредоточенной силой P (рису-
Рисунок 3 – Расчетная схема
жесткозащемленной по торцамнок 3). Были приняты следующие
цилиндрической панели, загруженной
исходныеданные:P  124,6 Н ;
сосредоточенной силой P
радиус цилиндраR  3,381 м ; толщина оболочки t  0,00476 м ; длина
образующих L  0,0254 м ; коэффициент Пуассона v  0,2 ; модуль упругости

E  7  10 4 МПа. Результаты расчётов оболочки в геометрически нелинейной
постановке представлены в таблице 2, в которой приводятся «физические»
значения нормальных напряжений  22 , а также величина прогиба в точке при-
Таблица 2 – Значения напряжений иложениясосредоточеннойсилы
прогиба при решении задачи в нелинейной
постановкеприразныхчислахшагов
нагружения. По данным таблицы 2,
можноутверждать,чтопри
увеличениичислашагов
нагружениянаблюдается
устойчивая сходимость вычислительного процесса, как по прогибу, так и по
напряжениям. В крайнем правом столбце таблицы 3 указано значение прогиба в
точке приложения сосредоточенной силы P , взятое из монографии Papenhausen
J. Различие в значениях прогиба не превышает 1%. Основываясь на анализе
табличных данных, можно сделать вывод, что разработанный алгоритм расчёта
позволяет получать приемлемые по точности значения параметров НДС тонкой
оболочки в геометрически нелинейной постановке с учетом сдвиговых
деформаций.
Пример 3. Был выполнен расчёт фрагмента оболочки, имеющей форму
эллиптическогоцилиндра,сжесткимзащемлениемпообразующей.
Интенсивность линейно распределенной вдоль образующей цилиндра нагрузки
q , приложенной на свободном краю, равнялась 1 кН / м (рисунок 4).
Параметры эллипса поперечного сечения имели следующие значения b  0,8 и
c  0,7 м. Толщина оболочки h  0,01 м, E  2  105 МПа, v  0,3 . Значение k
было принято равным k  0,5 . Число дискретных элементов равнялось 100.
Результаты расчётов сведены в таблицу 3, в которой приводятся значения
нормальных напряжений  22 на наружней  22н и внутренней  22в поверхностях
на свободном краю и в жесткой заделке, а также значения сдвиговых
напряжений  23 в сечении, находящемся вблизи жесткой заделки при
различных числах шагов нагружения. В крайнем правом столбце таблицы
приведенынапряжения,которыесоответствуютусловиямстатического
равновесия.
Таблица 3 – Численные значения напряжений  22
и  23 при действии линейно распределённой
нагрузки
СечениеНапря-Число шагов нагруженияКонтро
θ, раджение,406080100льное
МПазначени
е, МПа
Опор-в
 22-46.48-46.48-46.48-46.48-47.70
ноен
 2246.2646.2646.2646.2647.70
θ = 0,0в
 23-0.092-0.092-0.092-0.092-0.100
Рисунок 4 – Расчётная схемаСвобод-в
 220.0990.1010.1020.1020.100
ный
фрагмента эллиптическогокрай,
н
 220.0940.0960.0970.0970.100
цилиндра с жестким защемлениемθ = π/2
по образующей
По данным таблицы 3 можно отметить, что с увеличением числа шагов
нагружения наблюдается устойчивая сходимость вычислительного процесса.
Численные значения напряжений  22 и  23 на свободном краю при действии
линейно распределённой нагрузки q близки к значениям, полученным из
условия статического равновесия оболочки.
В четвертой главе представлен вывод новых интерполяционных
зависимостей для компонент вектора шагового перемещения и компонент
шагового вектора углов наклона нормали. Отличительной чертой этих
интерполяционныхзависимостейявляетсято,чтоинтерполяционные
зависимостибылискомпонованынедляотдельныхкомпонент,а
непосредственно для вектора шагового перемещения и шагового вектора угла
наклона нормали. После координатных преобразований были получены новые
интерполяционные зависимости, в которых компонента шагового вектора
перемещения и компонента шагового вектора углов наклона нормали
выражается через полный набор шаговых узловых неизвестных, а не отдельной
его части.
Интерполяционные выражения при векторной форме интерполяционной
процедуры записываются для самих векторов шагового перемещения и
шагового вектора углов наклона нормали
T 
w   wл    N wГ ;      ,
TT
(8)
11212111212121211441

гдеw   w w w w w
л Tijkli
,
  
…w,l w,i …w,l ;      ;
Tijkl

11241

w   w w w w w …w w
Г T

112
ijkli
,1
l
,1
i
,2

…w,l2 .

Шаговые векторы узловых точек конечного элемента и их производные
из (8) могут быть представлены компонентами узловых базисных векторов

w   u  a10   v  a20   w  a 0  ; w,  z11 a10   z12  a20   z1n a 0  ;1


w,  z 12 a10   z 22  a20   z 2n a 0  ;     1 a10    2  a20  ,
2(9)

где   i, j , k , l.
С учётом полученных выражений (9), соотношения (8) примут вид

w   N  Aw tw ;     A t .
TT

112
 
1212 1236 36114
 (10)
4881

где t    1i  2i  1 j  2 j  1k  2 k  1l  2l ,
T

18

t  u i v i wi u j v j w j …u l v l wl z11i z12i z1i z11 j z12 j z1j …z11l z 12i z 22i z 2i …z 2l .
T
w
136

   
Элементами матриц Aw и A в равенствах (10) будут соответственно
12 3648

 T   T 
подматрицы-строки aw0    a10  a20  a 0   и a0    a10  a20  , которые выражаются
1312

через орты декартовой системы координат

aw0    b  i ; a0    d  i .33312331
(11)
3121


С учётом i   b a 0  соотношения (11) могут быть представлены в виде
1 

313331
a   b b a ; a   d b a .
0
w


3333
10

31

0

23


33
10

31
(12)
3121

   
Матрицы Aw и A с учётом (11) можно представить в виде сумм
12 3648

A   a A  a A  a A ; A   a A  a A .
w
1236
w
1236
w
1236
w
1236

1
2
(13)
484848

Выполнив подстановки (11) и (12) в (10) соответственно, можно получить
равенства вида

ua10  va20  wa 0    a10 Aw1  a20 Aw2  a 0 Aw  Yw ;
T 
(14)
112 123612361236  361

 1a10   2 a20     a10 A1  a20 A2  Y ,
T (15)
14  4848  81

где Yw   G Tw wГ ; Y   T   1 2  .
3613636 36363618188
4141

Интерполяционные зависимости для компонент шагового вектора
перемещения и компонент шагового вектора угла наклона нормали можно
получить из (14) и (15)
u   Aw1 Yw ; v    Aw2 Yw ; w   Aw Yw ;
TTT

112123636111212363611121236361

 1    A1 Y ;  2    A2 Y .
TT
(16)
144881144881

Используяинтерполяционныевыражения(16),можнополучить
приращения деформаций на ( j  1) -м шаге нагружения в матричной форме

   BU ,

mn
544441
G
(17)
51

где  mn

  11 212 213  22 2 23 .
T

15

Пример 4. Была решена задача по определению НДС оболочки, имеющей
форму эллиптического цилиндра на которую действовала сосредоточенная сила
в середине пролёта, равнаяP . В эллиптическом сечении оболочки,
диаметрально противоположно действующей силе P (рисунок 5), находилась
шарнирная опора, которая препятствовала вертикальному смещению оболочки.
На торцах цилиндра была расположена скользящая заделка. Из-за наличия
плоскостей симметрии оболочки, расчет был выполнен для четвертой её части.
В качестве исходных данных были приняты следующие значения
E  0,738 105 МПа; v  0,3125; L  1,5 м. Отношение параметров эллипса
поперечного сечения эллиптического
цилиндра b / c было принято равным
1,5 / 1 .Результатыповариантных
расчётов представлены в таблице 4.
Анализируячисленныезначения
напряжений, приведённые в таблице 4
Рисунок 5 – Расчетная схема
эллиптического цилиндра,
можно убедиться, что различия между
загруженного сосредоточенной силой Pвариантами реализации интерполяции
шаговых узловых неизвестных весьма значительные. Исходя из условия
симметрии расчётной схемы, значения для нормальных и сдвиговых напряже-
Таблица 4 – Численные значения напряжений иний в точках A и B
прогиба эллиптического цилиндра, полученные при
решении задачи имеющей нелинейную постановку
должны совпадать, что и
Вариант интерполяционной процедурынаблюдаетсявовтором
МПа, прогиб,
Напряжения,
Координаты
точек x, θ

СтандартнаяВекторная
(м, рад)

вариантерасчета.
см

Число шагов нагружения

304050304050Полученныезначения
x = 0.0;
θ = 0.0в
5.075.075.073.773.783.78
напряженийвпервом
в
27.3627.3527.3525.3725.3625.36

 13в-1.46-1.48-1.49-1.04-1.06-1.07варианте расчета в точке A
 11н-4.89-4.87-4.86-3.54-3.53-3.52
больше напряжений, чем в
н
-31.99-31.99-31.99-30.25-30.26-30.26

 13н-1.46-1.48-1.49-1.04-1.06-1.07
точкеBвдвараза.
x = 0.0;
θ=π 11в2.892.892.893.773.783.78Значение прогиба в точке
 22в19.6819.6719.6725.3725.3625.36
приложения
в
-0.75-0.77-0.77-1.04-1.06-1.07

 11н-2.73-2.71-2.70-3.54-3.53-3.52
сосредоточенной силы P в
 22н-24.07-24.07-24.08-30.25-30.26-30.26первом варианте оказалось
н
-0.75-0.76-0.77-1.04-1.06-1.07
x = L/2
в 1,62 раз заниженным, чем
w-1.5252-1.5255-1.5256-2.0852-2.0858-2.0861
θ = 0.0
во втором варианте расчёта.
Анализ табличного материала позволяет сделать вывод, что при расчете
НДС эллиптических цилиндров в геометрически нелинейной постановке
общепринятая в МКЭ интерполяционная процедура приводит к некорректным
результатам.Прирасчётеэллиптическихцилиндроввгеометрически
нелинейной постановке следует применять разработанную векторную форму
интерполяции компонент вектора шагового перемещения и шагового вектора
углов поворота нормали. Следует отметить, что с увеличением кривизны
поперечного сечения эллиптического цилиндра отмеченная некорректность
стандартной интерполяционной процедуры будет нарастать. Выполнен расчёт
эллиптического цилиндра с вышеприведенными исходными данными при
соотношении параметров эллипса поперечного сечения равного 2 / 1 , а также
b / c  2,5 / 1 . Полученные результаты представлены в графической форме на
рисунках 6…9.

Рисунок 7 – Изменение значений
Рисунок 6 – Изменение значений
( 11в ) A( 22н ) A
11  в
вкоэффициентов  н
н
в зависимости
коэффициентовв зависимости от( 22 ) В
( 11 ) В
от отношения полуосей эллипса b/c
отношения полуосей эллипса b/c поперечного
поперечного сечения эллиптического
сечения эллиптического цилиндра
цилиндра

Рисунок 9 – Изменение значений
Рисунок 8 – Изменение значений
( )
в( 22н ) A
11в коэффициентов  н
н
в зависимости
( 22 ) В
11 A
коэффициентовв зависимости от22
( )
в
11 В
от отношения полуосей эллипса b/c
отношения полуосей эллипса b/c поперечного
поперечного сечения эллиптического
сечения эллиптического цилиндра
цилиндра
Введены коэффициенты, показывающие отношения нормальных и
сдвиговых напряжений в точках A и B при количестве шагов нагружения,
равным 50: 11в  ( 11в ) A /( 11в ) В ;  22н  ( 22н ) A /( 22н ) В . На рисунках 6…9 показаны

изменения вышеупомянутых коэффициентов 11в , 22н в зависимости от
соотношенияпараметров эллипсапоперечного сеченияэллиптического
цилиндра b / c в первом и во втором вариантах расчёта соответственно.
Как видно из рисунков, с увеличением отношения b / c значения
коэффициентов 11в ,  22н в первом варианте существенно возрастают и даже
изменяют свой знак, что является совершенно неприемлемым. Во втором
варианте наблюдаются стабильные значения коэффициентов 11в ,  22н , равные
единице, что и следует ожидать, исходя из условия симметрии расчётной схемы
эллиптического цилиндра.
На основании вышеизложенного табличного и графического материалов
можно сделать вывод о том, что при расчёте эллиптических цилиндров с
учетом сдвиговых деформаций в геометрически нелинейной постановке
необходимо использовать разработанную векторную форму интерполяционной
процедуры, применяемую как к компонентам вектора шагового перемещения,
так и к компонентам шагового вектора углов поворота нормали. Использование
стандартной интерполяционной процедуры в задачах исследования НДС
эллиптических цилиндров в геометрически нелинейной постановке приводит к
получению некорректных результатов.
ОСНОВНЫЕ РЕЗУЛЬТАТЫ И ВЫВОДЫ
1. Разработан алгоритм формирования матрицы жесткости и столбца
узловых усилий четырёхугольного конечного элемента на основе стандартной
для МКЭ интерполяционной процедуры и векторной формы интерполяции
искомых величин для анализа НДС тонких оболочек с учётом сдвиговых
деформаций в линейной постановке. На тестовых примерах доказано, что
векторная форма интерполяционной процедуры позволяет автоматически
учитывать смещения оболочечной конструкции как абсолютно твердого тела в
неявной форме.
2. Получены основные геометрические и физические соотношения на
шаге нагружения для тонкой оболочечной конструкции, определяющие связь
между приращениями деформаций и компонентами шагового вектора
перемещения, компонентами шагового вектора угла наклона нормали.
Сформированыматричныесоотношениямеждуконтравариантными
компонентамитензораприращенийнапряженийиковариантными
компонентами тензора приращений деформаций в произвольном слое оболочки
на( j  1) -ом шаге нагружения на основе предложенной гипотезы о
пропорциональностимеждукомпонентамидевиаторовприращений
деформаций и компонентами девиаторов приращений напряжений.
3.Разработаналгоритмформированияматрицыжесткости
четырёхузлового конечного элемента и столбца узловых усилий на шаге
нагруженияприреализацииразработаннойвекторнойформы
интерполяционной процедуры шаговых искомых величин.
4. Выполнен сравнительный анализ результатов конечно-элементных
решений оболочек с учетом сдвиговых деформаций при использовании двух
вариантов интерполяционной процедуры шаговых неизвестных, который
показал, что использование векторной формы интерполяции искомых величин
на шаге нагружения в конечно-элементном анализе оболочек с переменной
кривизной срединной поверхности позволяет получать корректные значения
прочностных параметров. Применение стандартной в МКЭ интерполяции
искомых величин при использовании криволинейных координат в расчетах
оболочек не приводит к корректным результатам.
5.Разработанныеалгоритмырасчетаоболочеквгеометрически
нелинейной постановке с учетом сдвиговых деформаций реализованы в пакетах
авторских программ по расчету на прочность оболочечных конструкций,
аппробированы и внедрены в практику инженерных расчетов.