Numerical modeling of nonlinear deformation processes for shells of medium thickness

Full Text

1. Введение Применение метода конечных элементов (МКЭ) к расчету тонкостенных конструкций в настоящее время является общепринятым и естественным. Разработке теоретических основ конечноэлементного анализа искривленных оболочек в линейной и нелинейной постановках посвящено значительное число научных публикаций, из которых отметим [1-5]. Среди большого разнообразия предложенных подходов для геометрически и физически нелинейных задач определенными преимуществами обладают так называемые двумерные и трехмерные изопараметрические конечные элементы (КЭ) тонкостенных конструкций [6-9]. Технология построения подобных элементов основана на использовании уравнений механики деформируемого твердого тела (МДТТ) и выполнении специфических гипотез в каждой квадратурной точке самостоятельно. Такой подход в полной мере позволяет применять апробированные методики нелинейного анализа поведения конструкций с учетом конечных деформаций [10-13], пластичности [14-16], накопления поврежденности материала [17], контактных задач [18-19], задач устойчивости [20-21], модели несжимаемых материалов [22] и т. д. В настоящей работе рассматриваются вопросы кинематики конечных деформаций и произвольных течений трехмерного континуума. Описываются тензоры деформаций и меры деформаций в разложении по основному и сопряженному базисам исходной и актуальной конфигураций. Построена математическая модель оболочечного конечного элемента путем введения в физическую модель материала Сетха. Рассмотрены числовые примеры. 2. Кинематика конечных деформаций Рассмотрим процесс деформирования оболочки в некоторой инерциальной системе отсчета, в которую введем декартовую систему координат с ортами e1, e2, e3. Пусть в исходной и актуальной конфигурациях радиус-векторы материальной точки имеют вид (1) Определим: - базисные векторы (2) - метрические тензоры (3) В дальнейшем нам понадобится вектор скорости v, который определим в следующем виде: (4) В задачах статики под вектором скорости понимается приращение радиус-вектора r по мере нарастания деформаций, то есть фактически v есть вектор приращения перемещений ΔU, где U - вектор перемещений. (5) Введем в рассмотрение группу тензоров, описывающих кинематику среды при ее деформации. Тензор градиента деформации (6) где Тензор градиента места (7) Обратный тензор градиента деформации (8) Обратный тензор градиента места (9) Отметим, что задание обратного тензора деформации в виде (8) может оказаться весьма неудобным, так как его компоненты вычисляются дифференцированием по координатам текущей конфигурации, которая часто является неизвестной. В этом случае целесообразно воспользоваться следующим представлением: (10) Обратный тензор градиента места определяется как транспонированный тензор. Далее при описании деформации использованы следующие соотношения: - правый тензор Коши - Грина (мера деформации Коши - Грина) (11) где - левый тензор Пиолы (мера деформации Альманси) (12) где - тензор деформации Коши - Грина (13) где - тензор деформации Альманси (14) где Следует отметить, что компоненты тензоров деформации Коши - Грина и Альманси в криволинейных базисах совпадают между собой. Справедливо соотношение, связывающее тензоры деформации Коши - Грина (13) и Альманси (14), (15) которое характерно для многих пространственных и материальных тензоров. Имеет место и обратное преобразование, то есть (16) Введены тензоры, используемые для описания течения среды. Базовыми тензорами здесь являются: - тензор пространственного градиента скорости (17) где (18) - тензор деформации скорости (19) где - тензор скорости поворота (20) где Материальная производная (полная производная по времени) тензора деформации Коши - Грина записана в виде (21) где Из (18) следует (22) Таким образом, скорость изменения тензора деформации Коши - Грина (21) и тензора деформаций скорости (19) в криволинейных базисах имеют одинаковые значения компонент. Следовательно, справедливо (23) Запишем (23) с учетом (15), (17) (24) Выражение в правой части (24) называют производной Ли (Lie Rate), для которой, с учетом (22), справедливо выражение (25) По аналогии с операцией дифференцирования по времени определим вариации основных тензоров: (26) 3. Физическая модель гиперупругого тела Для построения определяющих соотношений воспользуемся вторым уравнением термодинамики для изотермического деформирования упругой изотропной среды в виде (27) где ρ - плотность материала; ψ - функция свободной энергии; σ - тензор напряжений Коши; d - тензор деформации скорости, который определяется различными формами (19), (22), (23), (24), (25). Построим общие структурные соотношения, принимая в качестве базового выражения тензор деформаций Альманси (14). Выбор базового тензора деформаций предполагает, что функция свободной энергии зависит от компонент этого тензора. Таким образом, будем считать заданной функцию (28) По правилу дифференцирования скалярной функции по тензору имеем (29) где появляется тензор второго ранга Воспользуемся соотношением (24), из которого следует, что (30) Подставив (30) в (29) и далее в уравнение (27), получим (31) Преобразуем уравнение (31) к форме (32) Так как тензор скорости деформации d и тензор скорости вращения ω являются совершенно независимыми тензорами, то из (32) следуют два уравнения: (33) (34) Соотношения вида (33) называются физическими (определяющими) соотношениями, или уравнением состояния, и, по сути, определяют тензор напряжения в виде функции от тензора деформации. Уравнение (34) является ограничением на выбор функции свободной энергии (28). С учетом (34) определяющие соотношения (33) можно записать в виде (35) либо (36) Обе формы записи эквивалентны, вследствие справедливости (34). Для изотропного материала, свойства которого не зависят от направления, функция удельной потенциальной энергии деформации должна зависеть лишь от инвариантов соответствующих тензоров. Следовательно, выражение (28) упрощается до скалярной функции, зависящей от главных инвариантов (37) Здесь и далее обозначение скаляров упростим, а именно (38) В результате после несложных преобразований получим (39) где (40) Для получения представления тензора напряжений σ в виде суммы диад по соответствующему базису необходимо установить, в каком базисе определен тензор деформаций Альманси (14). Наиболее простой вид получаем при использовании базиса ei. В этом случае (41) и где образуется как матрица, обратная При использовании базиса текущей конфигурации получаем представление (42) и где - компоненты обратного тензора в базисе которые имеют достаточно сложный вид. В практической реализации вычисление компонент значительно проще, чем . Поэтому целесообразно вычислять через . Для этого из тождества (43) получим (44) Контравариантные компоненты тензора напряжений определяются в виде (45) Таким образом, представленные соотношения позволяют вычислять компоненты тензора напряжений в различных базисах при известных исходной и деформированной конфигурациях механической системы. 4. Разрешающее уравнение на шаге нагружения В качестве базового соотношения вариационное уравнение принципа виртуальных скоростей в текущей конфигурации можно записать в виде (46) где f* - вектор заданных внешних объемных сил; - вектор заданных напряжений на части поверхности Sσ, на которой определены силовые граничные условия. Технология вычислений представляет собой метод последовательных нагружений с определением текущей метрики как основной для вычислений. Итак, процесс деформирования представим как последовательность равновесных состояний Vk, которые реализуются при заданных значениях внешних сил . Определим в качестве основной неизвестной величины вектор скорости kυ, который можно трактовать как вектор приращения текущей конфигурации при переходе к новому состоянию Vk+1: (47) Разрешающее уравнение на текущем шаге строится путем линеаризации исходного уравнения (46) в предположении . В результате имеем (48) Чтобы полностью определить это уравнение необходимо построить выражение скорости напряжений для известной конфигурации kr через неизвестный вектор скорости (47) в виде линейной функции. Рассмотрим общий случай изотропного материала с определяющим уравнением (39). Индекс k для сокращения записи опустим. Справедливо (49) Далее распишем каждое из слагаемых: (50) Здесь использовались следующие представления: (51) Справедливо тождество (52) Если продифференцировать по времени тождество (52), то получим (53) Таким образом, удалось построить выражения всех слагаемых в соотношении (49) через два тензора скоростей Ȧ и ġ. Причем все тензоры, которые свертываются с ними, однозначно определяются текущей конфигурацией kr и могут быть вычислены по соответствующим формулам. Теперь необходимо выразить тензоры Ȧ, ġ через пространственный тензор градиента деформаций (17) в виде линейной зависимости. Из (24) имеем (54) то есть с учетом соотношения (19), (24) kȦ есть линейная функция, зависящая от kh. Из (14) получаем (55) Из (12) следует (56) Дифференцируя по времени тождество (57) получим (58) и подставим его в (56) (59) Собирая вместе представления (54), (55) и (59), получаем выражение (60) то есть соотношение, полностью аналогичное (54). Таким образом, представление для скорости изменения тензора напряжений (49) в глобальном базисе, то есть в виде (61) допускает представление (62) Выражение легко строится с помощью вышеприведенных соотношений. Следовательно, скорость изменения напряжений есть линейная функция от градиентов скоростей. Полученные соотношения представляют собой теоретическую основу конечноэлементного алгоритма исследования конечных деформаций нелинейно упругих тел при силовой нагрузке. Необходимо лишь добавить конкретную физическую модель в виде выражения функционала свободной энергии, справедливого для соответствующего материала. 5. Материал Сетха Расчет тонкостенных конструкций с учетом нелинейностей основывается на шаговых и итерационных методах. Выбор метода и алгоритма, реализующего его, зависит от типа нелинейности. В настоящей работе используется метод последовательных нагружений, который может быть естественно реализован в рамках МКЭ. Процесс деформирования представим в виде последовательности равновесных состояний V1, … Vk, Vk+1, … VN, где V1 и VN - области, занимаемые оболочкой в начальном и конечном деформированном состоянии, а Vk - произвольное промежуточное состояние. Действующая нагрузка достигается последовательным догружением на каждом шаге, причем количество шагов выбирается так, чтобы на каждом из них задача была квазилинейной. При такой постановке задача сводится к отысканию (k+1) - состояния при уже определенной геометрии и с накопленными напряжениями k-го состояния. Введем в рассмотрение текущую конфигурацию системы на k-м шаге нагружения и определим ее в следующем виде: (63) где (64) Соответственно вычисляем: - базисные векторы: (65) (66) - метрические тензоры (67) где (68) Если ввести в рассмотрение ковариантные компоненты метрических тензоров, то тензор деформации Альманси записывается следующим образом: (69) Так как компоненты тензоров деформаций Коши - Грина и Альманси в криволинейных базисах совпадают между собой, получаем тензор деформаций Коши - Грина (70) и тензор деформации Альманси (71) Введем в рассмотрение вектор приращения перемещений (72) где используются аппроксимации вида (64), то есть (73) Аналог тензора пространственного градиента скорости (74) будет представлен в виде (75) Симметричная часть этого тензора имеет вид (76) Аналогично можем записать вариации. Имеем (77) (78) Здесь имеют место следующие соотношения: (79) (80) Тензор истинных напряжений Коши определяется в виде (81) где введены ковариантные и контравариантные компоненты тензора напряжений. Запишем известное вариационное уравнение в скоростях напряжений (48) без учета массовых сил, предварительно сделав переход от и к приращениям и , приняв приращение времени Δt равным единице. Это уравнение будет иметь следующий вид: (82) В качестве физической модели используем материал Сетха, для которого справедлив закон Гука для тензора деформаций Альманси: (83) Распишем для приращения напряжений k-го состояния (84) где (85) После дискретизации (84) получаем (86) где (87) Далее по аналогии с [23] введем в (86) упрощенный закон Гука для приращений напряжений обжатия: (88) где E* - модуль жесткости на обжатие (в общем случае он может быть определен из экспериментальных данных). В частности, можно принять (89) Используя технологию усечения деформаций поперечного сдвига [23], получим (90) В целях уточнения физической модели, согласно теории, вводятся следующие соотношения: (91) где для принимаем соотношения (87), (90), а представим в виде (92) где, согласно технологии усечения деформаций поперечного сдвига, (93) что эквивалентно вычислению вращений в центре конечного элемента. Проведя некоторые преобразования, запишем матрицы геометрической жесткости второго и третьего слагаемых (82) в следующем дискретном виде: (94) В результате описанной выше конечноэлементной дискретизации получим систему линейных алгебраических уравнений (СЛАУ) (95) где - вектор приращения узловых перемещений; - матрица левых частей; - вектор приращения узловых сил; - вектор невязки. Решая систему линейных алгебраических уравнений (76), используя алгоритм переупорядочивания, известный как метод Катхилла - Макки [24], и определяя приращения перемещений, находим конфигурацию системы (l+1) и напряжения (96) 6. Числовые примеры Задача 1. Рассматривается тестовая задача изгиба полосы в кольцо. Исходя из кинематических соотношений вычислим и в узлах на свободном краю полосы с учетом Задача рассчитана с использованием предложенной выше методики. Длина полосы L = 200 см, толщина h = 1 см, ширина b = 5 см, модуль упругости Е = 20 000 кГ/см2, коэффициент Пуассона υ = 0. На рис. 1 изображено деформированное состояние полосы и несколько промежуточных этапов нагружения. Итог Рис. 1. Изгиб полосы в кольцо Figure 1. Bending the strip into a ring В случае, когда прилагаемая нагрузка разбивается на 1000 шагов нагружения, погрешность численного решения не более 1 %. Погрешность могла появиться за счет конечноэлементной аппроксимации исследуемой области. Приведенный числовой пример демонстрирует возможность настоящей методики в решении нелинейных задач теории оболочек. Задача 2. Рассматривается полусферическая оболочка с вырезом в полюсе под воздействием самоуравновешенной системы сил (растягивающих и сжимающих нагрузок в двух ортогональных радиальных направлениях). Из условий симметрии рассмотрим четверть полусферы (рис. 2). Изображение выглядит как рисунок, зарисовка, диаграмма Автоматически созданное описание Рис. 2. Полусферическая оболочка: радиус полусферы R = 10,0 см; толщина h = 0,04 см; модуль упругости Е = 6,825·107 кГ/см2; коэффициент Пуассона μ = 0,3; нагрузка F = 10λ кГ Figure 2. Hemispherical shell: radius of hemisphere R = 10.0 cm; thickness h = 0.04 cm; modulus of elasticity Е = 6.825·107 kg/cm2; Poisson's ratio μ = 0.3; load F = 10λ kg На рис. 3 изображен график максимальных перемещений Umax, см и Vmax, см при сетке 20×20 в сравнении с решениями других авторов [25-29]. На рис. 4 приведено деформированное состояние полусферической оболочки при λ = 16. Изображение выглядит как диаграмма, линия, График Автоматически созданное описание Рис. 3. График максимальных перемещений: Figure 3. Maximum displacement graph: Изображение выглядит как зарисовка, круг, рисунок, искусство Автоматически созданное описание Изображение выглядит как зарисовка, рисунок, диаграмма, искусство Автоматически созданное описание Рис. 4. Деформированное состояние полусферической оболочки Figure 4. Deformed state of a hemispherical shell Рис. 5. Цилиндрическая оболочка: радиус цилиндра R = 4,953 см; толщина h = 0,094 см; длина L = 10,35 см; модуль упругости Е = 10,5·106 кГ/см2; коэффициент Пуассона μ = 0,3125 Figure 5. Cylindrical shell: radius of cylinder R = 4.953 cm; thickness h = 0.094 cm; length L = 10.35 cm; modulus of elasticity Е = 10.5·106 kg/cm2; Poisson's ratio μ = 0.3125 Задача 3. Рассматривается растяжение цилиндрической оболочки путем приложения сосредоточенных сил P (рис. 5). Максимальное радиальное перемещение получено из аналитического соотношения На рис. 6 представлены радиальные перемещения на каждом шаге нагружения. Untitled1 Рис. 6. График радиальных перемещений Figure 6. Radial displacement graph Деформированное состояние цилиндрической оболочки при P = 1000 кГ представлено на рис. 7. Приближенные решения других авторов можно найти в [25-29]. цилиндр Рис. 7. Деформированное состояние цилиндрической оболочки Figure 7. Deformed state of a cylindrical shell 7. Заключение Представлена методика расчета напряженно-деформированных состояний тонкостенных конструкций сложной геометрии модифицированным трехмерным конечным элементом теории упругости с однослойной аппроксимацией по толщине и введен в рассмотрение материал Сетха. Приведенные численные расчеты тестовых задач демонстрируют работоспособность и достоверность предложенного алгоритма.

About the authors

Marat K. Sagdatullin

Kazan National Research Technological University

Author for correspondence.
Email: ssmarat@mail.ru
ORCID iD: 0000-0002-0535-4145

PhD in Physical and Mathematical Sciences, Associate Professor, Associate Professor of the Department of Basics of Design and Applied Mechanics

Kazan, Russian Federation

References

Supplementary files