Учет геометрической нелинейности в конечно-элементных прочностных расчетах тонкостенных конструкций типа оболочек

Полный текст

Введение 1 Современный анализ напряженно-деформи- рованного состояния (НДС) тонкостенных кон- струкций предполагает решение задачи в геометрически нелинейной постановке. При использовании численных методов рас- чета [1-6], в частности метода конечных элементов (МКЭ) [7-20], в решении нелинейных задач обычно используют шаговую процедуру нагружения [9; 10; 12; 14]. При этом возникает необходи- мость получения соотношений между приращениями деформаций, приращениями компонент вектора перемещения и приращениями их производных. В настоящей работе представлен вывод выше- упомянутых геометрических соотношений, включающих в себя деформации поперечного сдвига. Данные соотношения необходимы для формирования матрицы жесткости используемого конечного элемента на (j 1)-м шаге нагружения. 1. Геометрические соотношения При получении соотношений Коши на ( j 1)-м шаге нагружения последовательно рассматриваются три состояния оболочки: исходное и два деформированных - после j шагов нагружения и на (j 1)-м шаге нагружения. Исходное состояние описывается радиус-векторами Rr0 для точки 0 срединной поверхности и Rr0ζ для точки M0ζ , M находящейся на расстоянии ζ от срединной поверхности, причем r0ζ r0 r0 , (1) R  R ζen где ern0 - орт нормали к срединной поверхности в точке M0 . В процессе шагового нагружения точка M0ζ последовательно займет новые положения Mζ и M*ζ, определяемые соответствующими радиусвекторами: rζ r0ζ r r*ζ rζ r (2) R  R V; R  R W, где Vr и Wr - векторы перемещений точки M0ζ после j и ( j 1)-го шагов нагружения. При вычислении входящих в (2) векторов r и Wr можно воспользоваться одним из двух V вариантов. В первом варианте отсчет угла наклона нормали можно осуществить от ее исходного состояния [21]: V vr  r ζG W wr; r  r ζγ,r (3) где v v er ρ 0rρ vern0; w w er  ρ 0rρ wern0 - векторы перемещений точки M0 после j и ( j 1)-го шагов нагру- r ρ 0r ; γr  γρ 0erρ - векторы угла наклона жений; G  G eρ нормали после j и ( j 1)-го шагов нагружений ( ρ 1, 2 ). Во втором варианте отчет угла наклона нормали может осуществляться от ее деформированного состояния. В этом случае формулы (3) примут вид Vr   vr ζ ern Gr; Wr    wr ζ ern* γ ,r (4) где   ern ern ern0;   ern* ern* ern; ern и ern* - орты нормали после j и ( j 1)-го шагов нагружений. Ковариантные векторы базиса в трех рассматриваемых состояниях оболочки могут быть определены дифференцированием (1) и (2) по используемым глобальным криволинейным координатам. Например, если рассматривать в качестве рассчитываемой оболочки эллиптический цилиндр, то в качестве таких координат можно использовать осевую координату х и угловую координату θ: grα0  Rr,α0ζ; grα  Rr,αζ ; grα*  Rr,α*ζ, (5) где α последовательно принимает значения х и θ. Ковариантные компоненты тензора деформаций и тензора приращений деформаций после j шагов и на ( j 1)-м шаге нагружения могут быть получены из соотношений механики сплошной среды [22]: εζαβ gαβ  gαβ0  / 2;  εζαβ gαβ*  gαβ / 2. (6) Входящие в (6) ковариантные компоненты метрических тензоров в трех рассматриваемых состояниях могут быть определены скалярными произведениями (5): gαβ0  g g grα0  rβ0; αβ  g g gr rα  β; αβ*  g gr rα*  β*. (7) При использовании второго варианта отсчета угла наклона нормали (4) в соотношениях (6), выражающих приращения деформаций в произвольном слое оболочки через компоненты шагового вектора перемещения и компоненты шагового вектора угла наклона нормали, будут фигурировать как первые, так и вторые производные от компонент шагового вектора перемещения, что повлечет усложнение вычислительного алгоритма. Этим рассматриваемый вариант отличается от первого варианта отсчета угла наклона нормали, при использовании которого в соотношениях приращений деформаций фигурируют только первые производные от компонент вектора шагового перемещения. 2. Матрица жесткости на (j+1)-м шаге нагружения Элементом дискретизации выбирается четырехугольный фрагмент срединной поверхности с узлами в его вершинах. Столбцы шаговых узловых неизвестных в локальной 1 ξ, η 1 и глобальной х, θ системах координат будут иметь следующий вид: Wy Т w1yЛ Т w2yЛ      Т wЛy Т γ1 Т γ2 Т ; Л 1 44  1 12 1 12 1 12 1 4 1 4  WyГТ w1yГ Т w2yГ      Т wГy Т γ1 Т γ2 Т , (8) 1 44  1 12 1 12 1 12 1 4 1 4  где T qЛy  qq q qq q q q q q q qi j k l i,ξ ,ξj ,ξ ,ξ ,η ,η ,η ,ηk l i j k l ; 1 12 T qГy qq q qq q q q q q q qi j k l i,x ,jx ,kx ,lx ,θi ,θj ,θ ,θk l ; 1 12 T γρ γ γ γρi ρj ρkγρl; 1 4 q  w w w i j k l1, 2, ; , , , - узлы четырехугольного элемента дискретизации. Компонента шагового вектора перемещения и ее первые производные по глобальным координатам точки внутренней области конечного элемента аппроксимируются через узловые значения этой же компоненты с помощью интерполяционных выражений вида q  φ T qЛy ; q,α  φ,ξ T ξ,α  φ,η T η,α qЛy , (9) 1 12 12 1   где  φ T φ φ ...φ1 2 12 - матрица-строка, элементы 1 12 которой представляют собой произведение полиномов Эрмита третьей степени. Для компонент вектора угла наклона нормали γρ были использованы интерполяционные зависимости следующего вида: γρ  ψ Tγρy , (10) 1 4 4 1 где  ψ T ψ ψ ψ ψ1 2 3 4 - матрица-строка, элементы 1 4 которой представлены билинейными соотношениями локальных координат ξ, η . Функционал, выражающий равенство работ внешних и внутренних сил на ( j 1)-м шаге нагру- жения, записывается в виде T T П   εζαβ    σαβ  σαβdV      w  P  P dF , (11) V F T где εζαβ   ε11ζ 2 ε 12ζ 2 ε 13ζ εζ222 ε ζ23; 1 5 T σαβ   σ11σ12σ13σ22σ23 - приращение дефор- 1 5 маций и напряжений на ( j 1)-м шаге нагруже- T ния; σαβ σ σ σ σ11 12 13 22σ23 - напряжения, накоп- 1 5 ленные за j предыдущих шагов нагружения;  w T w w w1 2  - компоненты шагового вектора 1 3 перемещения точки срединной поверхности;  P T p p p1 2 3;  P T     p1 p2 p3 - внешняя по- 1 3 1 3 верхностная нагрузка за j шагов нагружения и ее приращения на ( j 1)-м шаге. Входящий в (11) столбец приращений контравариантных компонент тензора напряжений σαβ на основании закона Гука [22] может быть выражен через столбец приращений ковариантных компонент тензора деформаций εζαβ матричным способом σαβ С εζαβ, (12) 1 5 5 5 1 5 где С матрица упругости, при компоновке кото- 5 5 рой учтена общепринятая в теории оболочек [23] гипотеза о равенстве нулю нормальных напряжений, перпендикулярных срединной поверхности σ330. На основании соотношений (6) и аппроксимирующих выражений (9), (10) столбец приращений ковариантных компонент тензора деформаций εζαβ может быть выражен через столбец узловых неизвестных WyЛ в виде матричного произведения εζαβ  B WуЛ . (13) 5 1 5 44 44 1 С учетом (12), (13) и аппроксимирующих выражений (9), (10) функционал (11) примет следующий вид: П  WyГ TPRT     B T C B dV PR    WyГ  WyГ TPRT  B T σαβ dV  V V  WyГ TPRT  A T P dF  WyГ TPRT  A T P dF , (14) F F где матрица A определяется из равенства  w  A W yЛ. Выполняя над (14) процедуру минимизации T по WyГ , можно получить следующую систему алгебраических уравнений: PR     T  B T C B dV PR W  yГPR  T  A T  P dF V F  T T αβ T T  PR   B  σ dV PR   A  P dF , (15)  V F  которую можно записать в более компактном матричном виде: MG W yГ   f  NR , (16) где MG  PRT     B T C B dV PR  - матрица жестV кости конечного элемента на ( j 1)-м шаге нагружения;   f  PRT A T P dF - столбец узловых F усилий на ( j 1)-м шаге нагружения; NR - поправка Ньютона - Рафсона. 3. Пример расчета В качестве примера была решена задача по определению НДС цилиндрической панели, жестко защемленной по образующим и загруженной сосредоточенной силой P в середине пролета. При формировании матрицы жесткости и столбца узловых усилий конечного элемента использовались соотношения первого варианта отсчета углов наклона нормали как наиболее удобные с точки зрения организации вычислительной процедуры. Вследствие наличия плоскостей симметрии рассчитывалась 1/4 часть оболочки, которая представлялась одной полоской конечных элементов, ориентированных в кольцевом направлении. Были приняты следующие исходные данные: радиус цилиндра R  3,381 м; толщина оболочки t  0,00476 м; модуль упругости E  7 10 4МПа, коэффициент Пуассона v  0,2 ; величина сосредоточенной силы P 12,7 Н; длина образующих - 0,0254 м. Первоначально решалась задача в линейной постановке с целью установления необходимого числа элементов дискретизации. Таблица 1 Значения напряжений и прогиба при решении задачи в линейной постановке [Table 1. Values of stresses and deflection when solving a problem in a linear formulation] Координата Напряжение σ Количество элементов дискретизации [Number of sample elements] Известное решение точки, θ, рад МПа, прогиб ν, см [Known solution] [Point coordinate, θ, radian] [Voltage σ MPa, deflection ν, cm] 20 30 40 50 60 0,0, точка приложения силы P [0.0, point of application of force P] σ11в σ11н σв22 σн22 -1,67 2,68 41,17 -56,34 -1,30 2,31 43,55 -58,74 -1,06 2,07 44,82 -60,02 -0,91 1,92 45,63 -60,85 -0,8 1,81 46,23 -61,45 - - - - v -0,239 -0,240 -0,241 -0,241 -0,234 -0,241 0,128, жесткая заделка [0.128, hard fix] σ11в σ11н σв22 9,07 -11,78 45,30 9,94 -12,58 49,61 10,46 -13,09 52,21 10,83 -13,45 54,06 11,11 -13,74 55,48 - - - σн22 -58,99 -63,03 -65,55 -67,39 -68,82 - Таблица 2 Значения напряжений и прогиба при решении задачи в нелинейной постановке [Table 2. Values of stresses and deflection when solving a problem in a nonlinear formulation] Напряжение σ Число шагов нагружения [Number of sample elements] Известное решение МПа, прогиб ν, см [Voltage σ MPa, deflection ν, cm] 50 100 150 200 250 300 [Known solution] σв22 63,71 67,32 67,92 68,13 68,23 68,29 - σн22 -85,84 -84,96 -85,19 -85,34 -85,42 -85,48 - v -0,422 -0,429 -0,431 -0,432 -0,432 -0,432 -0,437 Результаты линейного расчета представлены в табл. 1: приведены численные значения нормальных напряжений на внутренней σв и наружной σн поверхностях оболочки в жесткой заделке, а также напряжения и прогиб в точке приложения сосредоточенной силы P в зависимости от числа элементов дискретизации. Как видно из табл. 1, при увеличении числа элементов дискретизации наблюдается сходимость вычислительного процесса как по напряжениям, так и по прогибу в точке приложения сосредоточенной нагрузки. В крайней правой колонке представлено известное линейное решение [24]. Из анализа данных табл. 1 можно сделать вывод, что разбиение 1/4 части исследуемой оболочки на 50 конечных элементов вполне достаточно, поэтому для решения задачи в нелинейной постановке было выбрано данное число элементов дискретизации. Результаты расчетов в геометрически нелинейной постановке представлены в табл. 2: приведены «физические» значения нормальных напря- жений σ22 и величина прогиба в точке приложения сосредоточенной силы в зависимости от числа шагов нагружения. Как видно из данной таблицы, с увеличением числа шагов нагружения наблюдается устойчивая сходимость вычислительного процесса, как по напряжениям, так и по прогибу. В крайней правой колонке приведено значение прогиба под сосредоточенной силой P, взятое из [24]. Вычисленное по разработанному в статье алгоритму значение прогиба оказалось зани- женным по сравнению с представленным в [24] всего на 1 %. Кроме того, следует отметить, что при увеличении количества элементов дискретизации и числа шагов нагружения величина прогиба будет монотонно возрастать. Заключение На основании анализа табличных данных можно сделать вывод, что разработанный алгоритм позволяет получать приемлемые по точности значения параметров напряженно-деформированного состояния тонких оболочек с учетом деформаций сдвига при расчете их в геометрически нелинейной постановке.

Об авторах

Юрий Васильевич Клочков

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

Автор, ответственный за переписку.
Email: klotchkov@bk.ru
SPIN-код: 9436-3693

д. т. н., профессор, заведующий кафедрой высшей математики

Российская Федерация, 400002, Волгоград, Университетский пр., 26

Анатолий Петрович Николаев

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

Email: klotchkov@bk.ru
SPIN-код: 2653-5484

д.т.н., профессор, профессор кафедры прикладной геодезии, природообустройства и водопользования

Российская Федерация, 400002, Волгоград, Университетский пр., 26

Тлек Рахметолович Ищанов

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

Email: klotchkov@bk.ru
SPIN-код: 1556-1368

к.т.н., старший преподаватель кафедры высшей математики

Российская Федерация, 400002, Волгоград, Университетский пр., 26

Александр Сергеевич Андреев

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

Email: klotchkov@bk.ru
SPIN-код: 7568-5011

старший преподаватель кафедры высшей математики

Российская Федерация, 400002, Волгоград, Университетский пр., 26

Михаил Юрьевич Клочков

Московский государственный университет имени М.В. Ломоносова

Email: klotchkov@bk.ru
SPIN-код: 2767-3955

студент 4-го курса физического факультета

Российская Федерация, 119991, Москва, Ленинские горы, 1

Список литературы