CC BY
46
Вычислительные технологии
Том 6, № 5, 2001
ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ C ОСОБЕННОСТЯМИ ДЕФОРМИРОВАНИЯ УПРУГОПЛАСТИЧЕСКИХ ОБОЛОЧЕК
ВРАЩЕНИЯ *
С. Н. КОРОБЕЙНИКОВ Институт гидродинамики им. М.А. Лаврентьева СО РАН
Новосибирск, Россия e-mail: [email protected]
The algorithm for numerical solution of quasistatic problems on axisymmetric deformation of shells of revolution, which are made from elastoplastic material, is suggested taking into account geometric nonlinearity. The determination of equilibrium configurations is reduced to the solution of Cauchy problem for a system of nonlinear ordinary differential equations due to finite-difference approximation of differential equations with respect to meridional coordinate. The special attention is given to determination of these configurations in neighborhoods of singular points of integral curve, namely, bifurcation and return points, where the matrix of system of equations degenerates. The auxiliary generalized problem on determination of eigenvalues and appropriate eigenvectors is suggested to be solved for continuation of the solution through singular points. The developed algorithm of solving problems in vicinities of singular points is applied to the problem on deformation of the longitudinally compressed simply supported circular cylindrical shell. The new solutions are obtained, which are verified by comparison with known experimental data. The singular points of both types (return and bifurcation) are obtained in the solutions of these problems. The singular points being simultaneously return and bifurcation points are also obtained.
Введение
В [1] приведен алгоритм численного решения нелинейных квазистатических задач об осе-симметричном деформировании оболочек вращения из упругого материала. Исходные нелинейные уравнения тонких оболочек Сандерса [2] дифференцируются по монотонно возрастающему параметру деформирования ¿. После конечно-разностной дискретизации уравнений по меридиональной координате получается система квазилинейных обыкновенных дифференциальных уравнений (ОДУ).
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, гранты №99-01-00525, 00-15-96180 и программы "Университеты России — фундаментальные исследования", грант №1795.
© С.Н. Коробейников, 2001.
Учет физической нелинейности, которая вводится через нелинейный закон связи скоростей компонент тензора напряжений со скоростями тензора деформаций (под скоростью понимается частная производная по £), приводит к тому, что результирующая система ОДУ становится нелинейной. В [3-5] предложены алгоритмы решения задач Коши для систем этих ОДУ методом пошагового интегрирования с итерационным уточнением решения.
Основные трудности, связанные с решениями задач Коши с нелинейными или квазилинейными системами ОДУ, возникают при достижении параметром £ некоторых критических значений, отвечающих вырождению матриц этих систем. Соответствующие точки интегральной кривой называются особыми. Общие алгоритмы численного решения нелинейных ОДУ, включающие особые точки интегральной кривой, представлены в [6, 7]. В настоящей работе предлагаются алгоритмы численного решения нелинейных ОДУ с особыми точками двух типов:
— бифуркации (ветвления) решений;
— поворота: в плоскости "характерное перемещение — параметр внешней силы" параметр внешней силы достигает локального максимума или минимума.
В [1, 8] также исследованы решения уравнений нелинейного деформирования в окрестностях особых точек другого типа. В этих точках достигаются предельные нагрузки: при стремлении параметра внешней силы к своему критическому значению (при котором вырождается матрица системы ОДУ) характерные перемещения обращаются в бесконечность. Предельные нагрузки достигаются также для тел из идеальных упругопластиче-ских материалов в решениях тех задач, в которых не учитывается геометрическая нелинейность [9]. При одновременном учете геометрической и физической нелинейностей (задачи такого типа рассматриваются в настоящей работе) предельные нагрузки, как правило, не достигаются, и особые точки этого типа здесь не рассматриваются.
Настоящая работа посвящена развитию алгоритмов и исследованию особенностей численных решений нелинейных уравнений, описывающих осесимметричное деформирование оболочек вращения из упругопластического материала в окрестностях особых точек, включающих как точки бифуркации, так и точки поворота. Для продолжения решения через точки бифуркации используется методика, развитая в [1] и адаптированная к настоящему классу задач. Для продолжения решения через точки поворота используется метод Бураго — Кукуджанова [10]1, который в [1, 8, 9] применен для определения равновесных конфигураций тел в окрестностях предельных состояний.
Предложенные алгоритмы апробируются на численных решениях задач об осесиммет-ричном деформировании упругопластических круговых цилиндрических оболочек, сжатых по оси. Решение тестовой задачи о потере устойчивости оболочки в упрощающем предположении о безмоментности докритического состояния указывает на надежность определения по разработанной вычислительной программе как равновесных конфигураций, так и форм потери устойчивости. В решениях этих задач с учетом моментности до-критического состояния на интегральных кривых решений получены особые точки обоих типов: в зависимости от геометрических параметров оболочки могут достигаться как точки бифуркации, так и точки поворота (локального максимума параметра внешней силы). Кроме того, получены качественно новые решения задач этого класса, в которых точки бифуркации совпадают с точками поворота. Некоторые из этих решений приведены в [8, 11].
хЭтот метод можно рассматривать как реализацию метода Рикса — Вемпнера [6], не нарушающую ленточность матриц ОДУ.
Опыт решения этих задач показывает необходимость тщательного исследования поведения решения в особых точках. Если решать задачу в "лоб", без выявления характера особых точек, то можно получить на первый взгляд парадоксальный результат: при вычислениях с простой и двойной точностями в одной и той же задаче по одному и тому же алгоритму можно получить качественно различные решения. Дело в том, что без решения вспомогательной задачи по определению собственных значений можно "пропустить" точки бифуркации, следуя по основной ветви решения. Также можно "пропустить" основную ветвь решения при выходе на боковую ветвь, решая задачу с простой точностью. Только исследование поведения решений в особых точках помогает выяснить то, что каждое из решений правильное, но они соответствуют разным ветвям, удовлетворяя одним и тем же уравнениям, граничным и начальным условиям задачи.
1. Уравнения квазистатического осесимметричного деформирования упругопластических оболочек вращения и алгоритм их решения
Уравнения, описывающие квазистатическое деформирование оболочек, формулируем относительно скоростей неизвестных величин. Под скоростями величин понимаем их частные производные по некоторому монотонно возрастающему параметру деформирования ¿, характеризующему процесс развития деформаций. Уравнения равновесия относительно скоростей получаются дифференцированием нелинейных уравнений Сандерса [2], описывающих поведение тонких оболочек в условиях малых деформаций, но, возможно, больших перемещений и умеренных поворотов. Определяющие соотношения материала оболочки задаются в виде однородных функций скоростей компонент тензора напряжений от скоростей компонент тензора деформаций. В силу того, что эти функции при нагружении и разгрузке материала оболочки записываются различными определяющими соотношениями, они относительно скоростей являются нелинейными. Используем теорию пластического течения с изотропным упрочнением материала (функция текучести Мизеса) [3-5, 12].
Уравнения, описывающие осесимметричное деформирование оболочек вращения, зависят только от одной пространственной переменной — меридиональной координаты. После конечно-разностной дискретизации уравнений по этой координате получаем задачу Коши вида [3-5]
А(Х, Х)Х = С, X = Х0 при г = 0. (1)
Здесь компонентами вектора X являются обезразмеренные меридиональная сила и1, перерезывающая сила с[ 1, изгибающий момент ш\, меридиональное перемещение и, нормальное перемещение т и угол поворота ф1 в узлах конечно-разностной сетки. Длина вектора неизвестных X равна N = 6М, где М — число узлов по меридиональной координате, включая граничные. В стандартной формулировке задачи вектор С предполагается известным (определяется скоростями изменения внешних сил). Матрица А имеет ленточную структуру, ее зависимость от текущих значений вектора Х обусловлена геометрической и физической нелинейностями уравнений, а зависимость от X является следствием нелинейности определяющих соотношений упругопластического материала. Точка над символом обозначает частную производную соответствующей величины по параметру г. Вектор X0 начальных значений неизвестных предполагается заданным.
Задача Коши (1) решается пошаговым интегрированием методом Эйлера с итерационным уточнением решения методом предиктор-корректор [3-5]. На каждой итерации система (1) решается относительно X методом матричной прогонки, следуя [13].
2. Особые точки интегральных кривых решений уравнений квазистатического деформирования
В дальнейшем предполагаем, что внешняя (консервативная) сила характеризуется одним скалярным параметром A(t) (интенсивностью) и некоторым постоянным вектором Co (ее распределением), так что вектор C можно представить в виде
C = A Co. (2)
В стандартной процедуре интегрирования системы (1) [3-5] параметр A предполагается заданной функцией A(t) (обычно A = t). Введем обозначение A(t) = A[X(t),X(t)]. Назовем особыми точками интегральной кривой такие, в которых
det A(t) = 0. (3)
Значения т параметра t, при которых выполняется равенство (3), назовем критическими.
Рассмотрим два возможных типа особых точек [6, 7]: точки поворота, в которых выполняется равенство
A = 0,
и точки бифуркации, в которых существует неединственное решение относительно X системы (1) (т. е. при t = т существует несколько продолжений решения). В последнем случае особая точка называется существенно особой [6, 7]. Допускаем также существование особых точек комбинированного типа, когда точка поворота одновременно является и точкой бифуркации.
В [1,8, 9] получены также решения задач с особыми точками, соответствующими предельным нагрузкам: таким критическим значениям параметра t, что
||X||2-► то,
t—>т
где ||X||2 = у X2 + + ... + X2 — евклидова норма вектора X. Значение Anm = lim A(t)
t—T
называется предельным. Отметим, что в решениях задач, рассматриваемых в настоящей работе, особые точки, соответствующие предельным нагрузкам, не встречаются.
При решении задач о деформировании оболочек из упругих материалов для обоих типов особых точек (бифуркации и поворота) достигаются собственные состояния [1], соответствующие таким равновесным конфигурациям, для которых существует нетривиальное решение однородной задачи
A[X(t)]W = 0 (4)
(при упругом деформировании оболочек матрица A от X не зависит). Определение собственных векторов W, составляющих базис нуль-пространства матрицы A, является важной частью алгоритма продолжения решения через точки бифуркации для получения боковых ветвей решений: решая задачу (1) с начальными данными, соответствующими
малым начальным неправильностям, близким по форме к собственному вектору W, получаем требуемую боковую ветвь [1]. Определение собственных состояний имеет также механический смысл: собственные состояния соответствуют смене устойчивых равновесных конфигураций на неустойчивые (по отношению к динамическим возмущениям) [12, 14]. Так как обычно основная ветвь решения послекритического режима деформирования соответствует неустойчивым равновесным конфигурациям, а боковая ветвь — устойчивым, больший практический смысл имеет определение боковой ветви. Кроме того, в точках бифуркации также теряется устойчивость квазистатических движений [12, 15]. Таким образом, для упругих оболочек потеря устойчивости квазистатических движений совпадает с потерей устойчивости равновесных конфигураций, сопровождаемой динамическим движением из неустойчивой равновесной конфигурации в устойчивую ("хлопок") [12].
Ситуация изменяется при решении задач о деформировании оболочек из упругопласти-ческого материала. В том случае, когда движение по боковой ветви решения происходит при возрастании параметра Л, задача о бифуркации решений не сводится к задаче об определении собственных состояний [16, 17]. При дискретизации уравнений, используемой в настоящей работе, задача (3) не сводится к задаче (4) вследствие зависимости матрицы А в (3) от X. Поэтому методику продолжения решения через особые точки, развитую в [1] и основанную на решении обобщенной задачи на собственные значения, для более общих классов задач, рассматриваемых в настоящей работе, прямо использовать нельзя.
На практике нас интересует наименьшее значение параметра Л, при котором происходит бифуркация решения. Для вычисления нагрузки, соответствующей этому значению Л (такая критическая нагрузка называется касательно-модульной [18]), надо в определяющих соотношениях упругопластичности при исследованиях поведения решения в точках бифуркации отбросить условие разгрузки по упругому закону и решить вспомогательную задачу по определению собственных состояний для "линейного тела сравнения" [16] (под этим телом понимается оболочка с материалом, полученным из упругопластическо-го материала при отбрасывании условий разгрузки). Критерий определения критической нагрузки, описанный выше, в [18] назван "критерием равноактивной бифуркации". При исследовании бифуркации переход от нелинейного тела к линейному телу сравнения на дискретном уровне означает игнорирование зависимости матрицы А от X, в этом случае задача (3) эквивалентна задаче (4).
Таким образом, для продолжения решения через точку бифуркации при достижении касательно-модульной нагрузки можно использовать собственные формы, полученные при решении задачи (4). Линеаризация задачи (4) рассмотрена в следующем разделе. Бифуркация решений задачи для оболочки из упругопластического материала в общем случае не сводится к решению задачи об определении собственных состояний. Это приводит к тому, что, в отличие от задач об упругом деформировании, потеря устойчивости квазистатического движения может происходить с устойчивыми равновесными конфигурациями [12, 18]. Впервые такой тип потери устойчивости наблюдал Ф. Шенли в экспериментах по выпучиванию стержней за пределом упругости [18]: после достижения критической нагрузки отклонение стержней, сжатых по оси, от прямолинейной формы происходит при продолжающемся нагружении без "хлопка".
3. Определение равновесных конфигураций оболочек в окрестностях особых точек
Если параметр Л предполагается заданной монотонно возрастающей функцией ¿, то равенство (3) в зависимости от вида матрицы А и вектора С соответствует тому, что либо решение системы (1) не существует (достигается точка поворота или предельная нагрузка), либо оно неединственно (происходит бифуркация решений задачи (1)) [6]. Эти варианты разрешимости или неразрешимости системы уравнений (1) следуют из альтернативы Фредгольма [19].
Решение задачи с помощью алгоритма, изложенного в разделе 1, с заданными значениями параметра Л(*) при достижении максимальной или предельной нагрузок получить нельзя. Тем не менее решить задачу по определению этих нагрузок и закритическому деформированию оболочек можно, модифицируя этот алгоритм с помощью метода Бураго — Кукуджанова [10]. В этом методе предлагается параметр Л ввести в число неизвестных, т. е. решать задачу не в Х-пространстве, а в (X, Л)-пространстве неизвестных. Основные положения этого метода применительно к рассматриваемому классу задач изложены в [1].
При численном решении задачи (1) (или модифицированной при использовании метода Бураго — Кукуджанова) критерием появления особых точек служит смена знака det А для двух соседних моментов времени ¿1 и ¿2 = ¿1 + Л* Исходя из критерия равноактивной бифуркации [18], вместо (3) рассматриваем задачу (4) (см. раздел 2). Линеаризуя задачу (4) на сегменте [¿1, ¿2], приходим к классической обобщенной задаче по определению собственных значений и соответствующих им собственных векторов [20]
(А* - ^В= 0, (5)
где матрицы А* и В* не зависят от параметра
Выражения для элементов матриц А* и В * получаются следующим образом. Представим матрицу А в виде суммы матриц
А(*) = А^) + А2[ф1(*),щ(*)], (6)
где элементы матрицы А1 зависят неявно от параметра * через текущие значения компонент тензора напряжений вследствие учета физической нелинейности, а матрица А2 появляется при учете геометрической нелинейности. Отметим, что для оболочки из упругого материала элементы матрицы А1 постоянны (не зависят от ¿) [1]. В рамках используемой теории оболочек Сандерса элементы матрицы А2 линейно зависят от ф1, п1. Тогда матрицу А2 можно представить в виде
А2 = ф1(£)А2 + П1(£)А2*, (7)
где А2, А2* — постоянные матрицы. Угол поворота ф1 и меридиональная сила п1 в общем случае зависят нелинейно от ¿. Предположим, что изменением элементов матрицы А1 на сегменте [£1,£2] можно пренебречь, так что
А1(£) = А^) V е М2]. (8)
Представим приближенно матрицу А в виде
д А
А(*) = А(*2) + т*(*)дА , (9)
где т* = Ь — Ь2 ^ 0 при Ь € [¿1, ¿2]. Из (6) - (9) получаем
дА дф дп , ,
— = А2 + А2*. (10)
дЬ дЬ 2 дЬ 2 v у
Представим приближенно
дф>1 _ Аф дП1 _ Ап (11)
дЬ _ АЬ ' дЬ _ АЬ ' ( )
где Аф1 = ф1(Ь2) — ф1(Ь1), Ап1 = п1(Ь2) — п1(Ь1). Пользуясь (9)-(11), аппроксимируем задачу (4) задачей (5), где
* .а _ у.
А* = А(*2)' В* = АФ1А2 + АщА2*' ^ = — А^ _ ^А^, 0 при Ь € (*1'*2). (12)
При геометрически линейном докритическом состоянии оболочки из упругого материала представления (8), (9), (11) являются точными, поэтому задачи (4) и (5) эквивалентны (линейная задача по определению критических состояний). При учете геометрической и/или физической нелинейностей задачу (5) называем линеаризованной задачей по определению критических состояний. Нас интересуют наименьшее собственное значение и соответствующий ему собственный вектор задачи (5). Они отыскиваются итерационным методом обратной степени [21].
После того как найдено собственное значение ^ст, приближенное критическое значение параметра Ь определяется из (12):
Ьсг _ ^2 — ^ст АЬ. (13)
Аналогично находится приближенное критическое значение Тст любой искомой величины Т (Ь):
Тст _ Т2 — ^ст(Т2 — Т1)' (14)
где Т1 _ Т(*1), Т2 _ Т(Ь2).
Численные расчеты, проведенные для оболочки из упругого материала, показывают [1], что собственные векторы и критические нагрузки, полученные при решении линеаризованной задачи (5) с учетом (14), определяются достаточно точно. Однако для оболочки из упругопластического материала допущение (8) довольно грубое. Поэтому для такого класса задач уточнять критическое значение Ьст € [Ь1'Ь2] надо не по формуле (13), а уменьшением шага интегрирования АЬ. В то же время собственные векторы задачи (4) находятся достаточно точно из решения линеаризованной задачи (5), что следует из решения тестовой задачи, приведенной ниже.
4. Тестовая задача по апробации алгоритма
определения равновесных конфигураций оболочек в окрестностях особых точек
Апробируем предложенный в разделе 3 алгоритм на задаче, имеющей аналитическое решение, — задаче об осесимметричном выпучивании продольно сжатой шарнирно опертой круговой цилиндрической оболочки в предположении о безмоментном докритическом состоянии. Ее решение сводится к определению наименьших корней нелинейного алгебраического уравнения, приведенного в [22].
<- L/2 -
Рис. 1. Шарнирно опертая круговая цилиндрическая оболочка, нагруженная осевой сжимающей силой N * > 0.
Рассмотрим круговую цилиндрическую оболочку (рис. 1) с радиусом срединной поверхности R, длиной L и толщиной h, сжатую по торцам осевой силой N * > 0 (силой, отнесенной к длине окружности радиуса R). Упрощающее (для получения аналитического решения) предположение о безмоментном докритическом состоянии вносит трудности в постановку задачи для получения численного решения, которое осуществлялось следующим образом. При решении задачи (1) об осесимметричном деформировании оболочки рассматривались граничные условия свободного от окружного стеснения края, воспроизводящие безмоментное состояние (т. е. ш\ = 0 при любом значении меридиональной координаты). Но при решении линеаризованной задачи (5) ставились граничные условия шарнирного опирания.
Материал оболочки предполагается упругопластическим с изотропным упрочнением с диаграммой одноосного деформирования вида [23, 24]
а/ау при а ^ ау, ,-1
еу \ n 1(а/ау)n + 1 — n 1 при а > ау,
где е, а — продольные деформация и напряжение соответственно; ау — начальное значение предела текучести; еу = ау/Е, Е — модуль Юнга; n — константа материала.
Расчеты проводились для оболочки со следующими геометрическими параметрами: R/h = 25, L/R = 2. Приняты механические константы материала: Е/ау = 250, n = 10, v = 0, 3 (v — коэффициент Пуассона). Задача (1) решалась с граничными условиями
n 1 = —n*, q1 = m 1 = 0 при £ = 0, (15)
i = ф 1 = ci1 = 0 при £ = 0.5,
где £ = s/L, s — меридиональная координата (рис. 1); n* = N*/(ауh); n1 = Ы11/(ауh) (N11 — погонное меридиональное усилие в оболочке); u = U/h (U — меридиональное перемещение). Краевые условия (15) при £ = 0 соответствуют свободному относительно
нормальных перемещений краю, подвергающемуся осевой сжимающей силе, а при £ = 0.5 — условиям симметрии решения относительно плоскости, проходящей через середину образующей оболочки перпендикулярно к ней. По меридиональной координате £ сегмент [0; 0.5] разбивается на 25 равных интервалов. За параметр деформирования принимается смещение левого торца оболочки, т. е. £ = и |^=0. Интегрирование системы (1) проводилось с шагом Д£ = 0.005. Начальные условия в (1) принимались нулевыми (т.е. рассматривается оболочка без начальных неправильностей). Численное решение задачи (1) доводится до £тах = 0.385.
Задача (1) (основная задача) решается для моделирования безмоментного докритиче-ского состояния оболочки. Для определения критических нагрузок и форм потери устойчивости шарнирно опертой оболочки рассматриваются однородные граничные условия, полученные из соответствующих условий шарнирного опирания (см. раздел 5):
п 1 = го = т 1 = 0 при £ = 0, „Л
п I (16)
и = го = т 1 = 0 при £ = 1,
где г = Ж/Л,, Ш — нормальное перемещение (прогиб). Отметим, что формы потери устойчивости могут быть как симметричными, так и антисимметричными относительно середины образующей оболочки, поэтому здесь рассматривается не половина, а весь сегмент [0; 1] координаты £, который разбивается на 50 равных интервалов. При решении задачи (5) компоненты тензора напряжений и распределение п1 и ф1 при £ £ [0; 0.5] брались из решения основной задачи, а при £ £ (0.5; 1] восстанавливались из условий симметрии.
При численном решении основной задачи (1) с граничными условиями (15) получено, что А вспомогательной задачи с граничными условиями (16) на сегменте £ £ [0; 0.385] два раза меняет знак: на подсегментах [£1 = 0.345; £2 = 0.350] и [£3 = 0.370; £4 = 0.375]. Из решения основной задачи получены следующие значения осевой сжимающей силы п* на границах этих подсегментов:
п*(£1) = 1.382; п*(*2) = 1.385; п*(£з) = 1.395; п*(£4) = 1.398.
Таким образом, два нижних критических значения осевой сжимающей силы п^1 и п^2, найденные при численном решении рассматриваемой задачи, заключены в интервалах
пП1 £ (1.382; 1.385), <2 £ (1.395; 1.398). (17)
Обозначим полученные в аналитическом решении [22] критические значения осевой силы п* и собственные формы соответственно через п^ и
г
^т кп£ (г = 1, 2,...), (18)
где кг — положительные целые числа. В зависимости от того, являются ли кг четными или нечетными, собственные функции г становятся соответственно антисимметричными или симметричными относительно середины образующей оболочки. При осесимметричном выпучивании, описываемом уравнениями Сандерса, получаем [22] набор трех нижних значений критической силы п* и соответствующих им собственных форм (18):
<1 = 1.383, <2 = 1.393, паз = 1.410, ( )
к1 = 5, к2 = 6, кз = 4. ( )
Сравнивая (17) и (19), отмечаем, что границы интервалов, в которых заключены п^1 и <2, определены достаточно точно. При первой смене знака А в интервале (£1,£2)
собственная функция weifl)1, построенная из компонент собственного вектора W (полученного из решения задачи (5)), с высокой точностью совпадает с функцией w1 в (18) при k1 = 5. Использование формулы (14) для уточнения критической силы дает значение n*n1 = 1.372. При второй смене знака det A в интервале (t3,t4) аналогично получаем n*n2 = 1.387, а собственная функция weig,2 близка к w2 с k2 = 6. Отсюда следует, что решение задачи (5) помогает определить только собственные векторы Wi, представляющие базис нуль-пространства матрицы A при достижении особых точек решения уравнения (1). Критические значения параметра деформирования и соответствующие им критические значения внешних сил уточнить не удается (при "уточнении" этих значений по формулам (14) значения n^1 и n*n2 даже выходят за границы интервалов, в которых они должны находиться). Это заключение, сделанное на основе решения тестовой задачи, согласуется с заключением, сделанным в [20].
5. Численные решения задач о деформировании продольно сжатых шарнирно опертых круговых цилиндрических оболочек из упругопластического материала
5.1. Постановка задачи
Определим закономерности деформирования и потери устойчивости оболочки из упру-гопластического материала, сжатой по оси, в зависимости от изменения геометрических параметров R/h и L/R по алгоритму решения нелинейных задач теории оболочек, изложенному в разделе 3 и включающему определение равновесных конфигураций в окрестностях особых точек. Результаты расчетов представлены на рис. 2-6 в осях "||w|| ~ n*", где
На этих рисунках также приведены равновесные конфигурации в осях " £ ~ ад" в некоторых характерных точках интегральных кривых. Около каждой из конфигураций дано распределение зон активных пластических деформаций по сечению оболочки (закрашены). Также на каждом рисунке показаны собственные формы (£), полученные при решении задачи (5). Около этих форм приведены распределения зон пластических деформаций с игнорированием условий разгрузки по упругому закону (см. раздел 2).
Предполагается, что левый торец оболочки подвергается осевой сжимающей силе, а правый закреплен относительно осевых перемещений (см. рис. 1). На левом торце оболочки ставится граничное условие шарнирного опирания
либо на середине образующей (условие симметрии деформирования относительно плоскости, проходящей через середину образующей перпендикулярно к ней)
п 1 = —п*, ад = т 1 = 0 при £ = 0. Второе граничное условие ставится либо на правом торце (шарнирное опирание)
й = ад = т 1 = 0 при £ = 1,
(20)
(21)
U = ф 1 = = 0 при £ = 0.5.
(22)
Задаются следующие начальные условия интегрирования системы (1)2:
(23)
где амплитуда 5 и целое число полуволн m задаются, b = h/L. Значения 5 = 0 и m = 0 соответствуют оболочке с начальным несовершенством (начальной неправильностью), а 5 = m = 0 — идеальной оболочке.
Материал всех рассматриваемых далее оболочек предполагается упругопластическим с теми же характеристиками, что и в разделе 4. В качестве параметра деформирования t для всех оболочек использовано значение прогиба в районе действия краевого эффекта. При бифуркации решений ветвь, соответствующую симметричному деформированию относительно середины образующей, назовем основной ветвью решения, а ветвь, на которой эта симметрия теряется, — боковой ветвью. Расчеты проводились как с простой, так и двойной точностью. Мантисса вещественных чисел в первом случае составляет шесть десятичных разрядов, а во втором — шестнадцать.
Представим полученные решения.
5.2. Решение для оболочки с R/h = 25, L/R = 0.5
Численные решения для оболочки с этими параметрами приведены на рис. 2. Здесь t = w 1^=0.26, At = 0.003. В этой задаче как для идеальной оболочки, так и для оболочек с начальными неправильностями получена одна особая точка, являющаяся точкой поворота (в этой точке достигается максимальная нагрузка). Как видно из графиков, приведенных на рис. 2, влияние начальных неправильностей антисимметричной формы (m = 2) на решение задачи для идеальной оболочки оказывается несущественным.
5.3. Решение для оболочки с R/h = 25, L/R = 0.7
Решения для этой оболочки приведены на рис. 3. Так же, как и в предыдущем подразделе, t = w |g=0.26, At = 0.003. В решении для идеальной оболочки обнаружена точка бифуркации. Основная ветвь решения (ветвь 1 на рис. 3) получена с граничными условиями как (20), (21), так и (20), (22). В первом случае решение найдено при вычислениях с двойной точностью. Решение для боковой ветви (ветвь 2) получено как при решениях задач для оболочек с начальными неправильностями малой амплитуды с m =2 и 5 = 10-12, 10-8, так и при решении задачи (1) для идеальной оболочки при вычислениях с простой точностью. Последний расчет можно рассматривать как решение задачи с начальной неправильностью малой амплитуды и формой случайного вида, вводимой погрешностью вычислений на компьютере.
При расчетах, проведенных с двойной точностью для идеальной оболочки, в решении получены последовательно две особые точки: сначала точка бифуркации, а потом точка поворота, соответствующая максимальной нагрузке. На интегральной кривой, соответствующей боковой ветви, получена только одна особая точка — точка поворота, соответствующая максимальной нагрузке.
2Точнее говоря, с помощью (23) получаются компоненты вектора Xq в (1).
Рис. 2. Поведение решения задачи о деформировании оболочки с геометрическими параметрами
R/h = 25, L/R = 0.5: А — начало пластического течения; х — точка максимума;--решения
для идеальной (простая и двойная точности, условия симметрии) и неидеальной (wo = ö sin2n£; ö = 10-8, 10-4) оболочек;---решение для неидеальной оболочки с wo = 0.01sin2n£.
Точка бифуркации соответствует потере устойчивости квазистатических движений (см. раздел 2). Отсюда и следует чувствительность решения к начальным несовершенствам, имеющим форму собственных функций [12, 15]3. Симметричная собственная форма, соответствующая точке поворота на основной ветви решения, на рис. 3 не приводится. Дана только антисимметричная собственная форма, полученная в точке бифуркации. Переход решения при движении вдоль интегральной кривой на боковую ветвь сопровождается локальной областью разгрузки материала оболочки, которая возникает в послекритическом режиме деформирования сразу же за точкой бифуркации. Эта область постепенно расширяется, в некоторый момент деформирования она приведена на рис. 3. Такой характер смены области пластического нагружения на упругую разгрузку качественно согласуется с аналогичным поведением зон нагрузки и разгрузки в решении задачи о выпучивании стойки из упругопластического материала [24, 25].
3В расчетах использовались начальные неправильности, согласующиеся с собственными формами в районах действия краевых эффектов.
5.4. Решение для оболочки с R/h = 25, L/R = 2
Графики расчетов при t = w |^=0.08, At = 0.0025 представлены на рис. 4. Продолжение решения за точкой бифуркации вдоль основной ветви (ветвь 1 на рис. 3) получено как при расчетах с граничными условиями (20), (22) с выполнением условий симметрии, так и с граничными условиями (20), (21) при вычислениях с двойной точностью для идеальной оболочки и для оболочек с начальными неправильностями с m =2 и 8 = 10-12, 10-8, 10-6, 10-5.
При получении решения для идеальной оболочки критерий определения особых точек по смене знака det A "сработал" только в решении задачи с граничными условиями (20), (22). При решении задачи с граничными условиями (20), (21) смены знака det A не произошло. Причиной послужило то обстоятельство, что при решении задачи (5) с граничными условиями (16) нижнее собственное значение имеет кратность 2. Поэтому в этой и последующих далее задачах (см. подразделы 5.5, 5.6) задача (5) решается с однородными
Рис. 3. Поведение решения задачи о деформировании оболочки с геометрическими параметрами К/Н = 25, Ь/К = 0.7: Д — начало пластического течения; х — точка максимума; □ — точка бифуркации решений; ветвь 1 — решение с двойной точностью для идеальной оболочки; ветвь 2 — решение с простой точностью для идеальной оболочки и с двойной точностью для неидеальной
оболочки с Шо = 5 8ш2п£, 5 = 10-12, 10-8;---решение для неидеальной оболочки с Шо =
10-4 8ш2п£;---решение для неидеальной оболочки с Шо = 0.01 8ш2п£.
граничными условиями
и 1 = т = т 1 = 0 при £ = 0, (24)
а при £ = 0.5 выполняются либо условия симметрии
и = ф 1 = & = 0, (25)
либо антисимметрии
и = т = т 1 = 0. (26)
Когда в решении задачи с граничными условиями (20), (22) меняет знак А, решались две задачи (5). При решении первой задачи использовались граничные условия (24), (25), а при решении второй — (24), (26). В решениях обеих задач получено одно и то же наименьшее собственное значение ^сг, что эквивалентно наименьшему собственному значению кратности 2 при решении задачи (5) с граничными условиями (16). Симметричная т3 и антисимметричная та собственные формы для полной оболочки восстанавливаются из условий симметрии и антисимметрии (см. рис. 4). Таким образом, нуль-пространство матрицы А, полученной в особой точке интегральной кривой решения задачи (1) с граничными условиями (20), (21), порождается двумя линейно независимыми векторами W1 и (являющимися собственными векторами задачи (5)), из компонент которых определяются собственные формы т3 и та.
Наличие кратного собственного значения и вид собственных векторов указывают на то, что в настоящей задаче реализуется ситуация, на возможность которой указано в [18, 26], когда максимальная нагрузка достигается одновременно с гладкой бифуркацией (бифуркацией порядка выше первого4). Задачи этого же класса, когда бифуркация решений происходит в точке поворота, приведены в [23] для разного типа конструкций из упругопластического материала (арки, упругого стержня на упругопластическом основании, сжатой пластинки). Такого же типа решения получены в задачах о деформировании и потере устойчивости атомных решеток [27]. Некоторые задачи, в которых встречается бифуркация решений порядка выше первого, рассмотрены в [28].
Из анализа графиков, приведенных на рис. 4, следует, что оболочка с Ь/Я =2 менее чувствительна к начальным несовершенствам, чем оболочка с Ь/Я = 0.7 (см. рис. 3). Решение следует вдоль боковой ветви как при расчетах с двойной точностью, так и при расчетах с простой точностью для оболочек с начальными неправильностями с т = 2 и довольно большими амплитудами5 8 = 10-3, 10-4.
Из вида представленных на рис. 4 равновесных конфигураций оболочек следует, что ветвь основного решения в послекритическом режиме деформирования характеризуется ростом кольцеобразных складок около обоих торцов, сопровождающихся активным пластическим деформированием. На боковой ветви решения такая складка продолжает расти только около одного из торцов. Разгрузка материала около другого торца происходит в конечной области сразу же при выходе решения на эту ветвь, а не постепенно, как для оболочки с Ь/Я = 0.7 (см. подраздел 5.3), для которой происходит бифуркация решения при продолжающемся нагружении. Такая же смена конечной области активного
4Касательные к интегральной кривой в точке бифуркации для основной и боковой ветвей решения совпадают [18].
5Напомним, что решения для оболочек с начальными неправильностями меньшей амплитуды следуют вдоль основной ветви.
О 0.05 0.1 |М
Рис. 4. Поведение решения задачи о деформировании оболочки с геометрическими параметрами К/Н = 25, Ь/К = 2: Д — начало пластического течения; 1X1 — точка максимума, совпадающая с точкой бифуркации решений; ветвь 1 — решение с двойной точностью для идеальной оболочки и неидеальной оболочки с -—о = 5 8ш2п£, 5 = 10-12, 10-8, 10-6, 10-5; ветвь 2 — решение с простой
точностью для идеальной оболочки и неидеальной оболочки с шо = 10-4 8ш2п£;---решение
для неидеальной оболочки с -о = 10-3 8ш2п£;---решение для неидеальной оболочки с
ш0 = 0.01 8Ш2П{.
пластического деформирования на соответствующую область с упругой разгрузкой материала получена в задаче о деформировании стойки за пределом упругости при достижении приведенно-модульной нагрузки [25]6. В настоящей задаче бифуркационная нагрузка является приведенно-модульной, так как достигается в точке поворота при п* = 0. Эта нагрузка одновременно является и касательно-модульной в силу того, что при ее достижении впервые происходит бифуркация решений.
5.5. Решение для оболочки с И/Н = 50, Ь/И = 2
Графики расчетов при £ = т |^=0.0б, = 0.005 даны на рис. 5. Качественное поведение решения остается таким же, как и в подразделе 5.4, но кольцеобразные складки локали-
6 Приведенно-модульной нагрузкой [18] называется нагрузка, соответствующая точке бифуркации, для которой боковая ветвь решения в начальном послекритическом режиме деформирования характеризуется постоянной внешней силой (А = 0), так что вектор С = 0 в (1) в соответствии с (2).
зованы в большей степени около торцов оболочки.
5.6. Решение для оболочки с R/h = 100, L/R = 2
Численные решения отображены на рис. 6. Принимается t = w |j=0.04, At = 0.01. Качественное поведение решения остается таким же, как и в решениях, приведенных в п. 5.4, 5.5, с еще большей локализацией кольцеобразных складок. Но в отличие от поведения решений, приведенных в предыдущих подразделах, здесь при вычислениях с простой точностью в послекритическом режиме деформирования решение следует вдоль основной ветви.
Исследовано влияние формы начальных неправильностей на поведение решения. Рассматривались антисимметричные начальные неправильности с m = 2 и m = 12. Из приведенных на рис. 6 графиков следует, что для одинакового значения амплитуды 5 большее влияние оказывают начальные неправильности с m = 12, соответствующие форме потери устойчивости с безмоментным докритическим состоянием (см. [22] и раздел 4 настоящей статьи). Эти начальные неправильности больше согласуются в районе действия краевого эффекта с формой weig = wa (рис. 6) и вносят большие возмущения в решение. В [1] показано, что при деформировании упругой оболочки наибольшее влияние на поведение решения в окрестностях особых точек оказывают начальные неправильности, совпадающие по форме с собственными функциями. Из настоящих расчетов следует, что аналогичный эффект наблюдается и в решении для оболочки из упругопластического материала.
5.7. Обсуждение результатов расчетов
В теоретических и экспериментальных исследованиях установлено [29, 30], что выпучивание круговых цилиндрических оболочек, сжатых по оси, при упругом деформировании происходит по неосесимметричной форме. Поэтому решения задач о деформировании и потере устойчивости оболочек из упругого материала по осесимметричной форме (не реализуемой на практике), приведенные в [1], имеют главной целью отработку алгоритма по решению задач деформирования оболочек из упругого материала в окрестностях особых точек.
Положение изменяется при решении подобного класса задач о деформировании оболочек из упругопластического материала. Для круговых цилиндрических оболочек, сжатых по оси, при деформировании за пределом упругости реализуется осесимметричная форма выпучивания [22, 29, 30] для оболочек средней длины7, являющихся предметом настоящего исследования. Критические значения осевой сжимающей силы, полученные в предположении о безмоментном докритическом состоянии [22]8, достаточно близки к соответствующим значениям силы, полученным в настоящей работе (см. рис. 4-6). В то же время сильно различаются формы выпучивания: синусоида — в решении с допущением безмоментности докритического состояния [22]; локальное выпучивание около одного или обоих торцов — в настоящем решении.
В экспериментах по сжатию круговых цилиндрических оболочек с геометрическими параметрами, близкими к рассмотренным в подразделах 5.4-5.6, происходит потеря устой-
7При достаточно большом значении параметра L/R потеря устойчивости оболочки происходит по стержнеобразной (неосесимметричной) форме [22].
8Эти значения (обозначенные через n*mb) отмечены горизонтальными штриховыми линиями на рис. 4-6.
п
«1=1.220
1.0
0.75 -
0.5 -
0.25 -
0
Рис. 5. Поведение решения задачи о деформировании оболочки с геометрическими параметрами R/h = 50, L/R = 2: Д — начало пластического течения; И — точка максимума, совпадающая с точкой бифуркации решений; ветвь 1 — решение с двойной точностью для идеальной оболочки и для неидеальной оболочки с wo = 10-8 sin2n£; ветвь 2 — решение с простой точностью для идеальной оболочки и с двойной точностью для неидеальной оболочки с wo = S sin 2п£ с амплитудами S = 10-6, 10-4;---решение для неидеальной оболочки с wo = 0.01sin2n£.
чивости за пределом упругости по осесимметричной форме. При этом реализуется как несимметричная относительно середины образующей форма выпучивания (см. рис. 14.1 в [29], рис. 1 в [31], рис. 11 в [32], рис. 6 в [33], рис. 4.30 в [34]), так и симметричная (см. рис. 4.30, 4.32 в [34]). Формы потери устойчивости, полученные в экспериментах и приведенные в упоминаемых выше работах, так же, как и полученные в настоящем исследовании (см. подразделы 5.4-5.6), — кольцеобразные складки около одного или обоих торцов.
Рис. 6. Поведение решения задачи о деформировании оболочки с геометрическими параметрами R/h = 100, L/R = 2: А — начало пластического течения; И — точка максимума, совпадающая с точкой бифуркации решений; ветвь 1 — решения с двойной и простой точностями для идеальной оболочки и неидеальной оболочки с wo = ö sin2n£, ö = 10-8, 10-4, 10-3; wo = 10-6 sin 12n£; ветвь 2 — решения с двойной точностью для неидеальных оболочек с начальными неправильностями вида w0 = ösin2n{, ö = 0.005, 0.01; w0 = ö sin12n{, ö = 10-5, 10-4;---решение
для неидеальной оболочки с wo = 0.001 sin 12п£;---решение для неидеальной оболочки с
wo = 0.01 sin 12п£.
Заключение
Развитый в работе алгоритм позволяет определить решения рассматриваемого класса задач в окрестностях особых точек интегральной кривой. Решение вспомогательной задачи на собственные значения позволяет найти боковые ветви решения основной задачи, повторяя численное решение с начальными неправильностями, учитывающими вид найденных собственных векторов.
Полученные решения задач по потере устойчивости круговой цилиндрической оболочки из упругопластического материала показывают, что при вычислениях с простой точностью точку бифуркации можно вообще пропустить. В этом случае получаем решение
сразу для боковой ветви. Наоборот, решая задачу с двойной точностью (без рассмотрения возможности появления особых точек на интегральной кривой), можно получить решение только для основной ветви. Таким образом, задачи с особыми точками интегральной кривой требуют тщательного подхода как к разработкам алгоритмов их решения, так и к самим решениям. Кроме того, требуется проверка устойчивости этих решений.
Для круговой цилиндрической оболочки из упругопластического материала, сжатой по оси и шарнирно опертой по торцам, в зависимости от геометрических параметров в решениях достигаются особые точки интегральной кривой следующих типов:
— точка поворота, соответствующая потере устойчивости равновесных конфигураций;
— точка бифуркации с продолжающимся нагружением на боковой ветви решения. В этой точке происходит потеря устойчивости квазистатических движений, но равновесные конфигурации остаются устойчивыми. При дальнейшем нагружении достигаются точки поворота на каждой из ветвей с потерей устойчивости равновесных конфигураций;
— точка поворота с одновременной бифуркацией решений. Теряется устойчивость квазистатического движения с одновременной потерей устойчивости равновесных конфигураций.
Список литературы
[1] Коробейников С. Н. Численное решение нелинейных задач о деформировании упругих оболочек вращения в собственных состояниях // Сиб. журн. вычисл. математики. 2000. Т. 3, №1. С. 43-56.
[2] Sanders J. L. Nonlinear theories for thin shells // Q. Appl. Math. 1963. Vol. 21, No. 1. P. 21-36.
[3] Волчков Ю.М., Коробейников С.Н. Численное решение упругопластических задач теории оболочек // Численные методы решения задач теории упругости и пластичности: Мат. 5-й Всесоюз. конф. Ч. 2 / АН СССР. Сиб. отд-ние. Ин-т теоретической и прикладной механики. Новосибирск, 1978. С. 40-47.
[4] Волчков Ю.М., Коробейников С.Н. Численное определение предельных и бифуркационных нагрузок упругопластических оболочек вращения // Динамика сплошной среды / АН СССР. Сиб. отд-ние. Ин-т гидродинамики. Новосибирск, 1980. Вып. 45. С. 58-78.
[5] Волчков Ю.М., Коробейников С.Н. Оценка предельной нагрузки упругопластических оболочек вращения // ПМТФ. 1981. №4. С. 146-150.
[6] ШАллшилин В. И., Кузнецов е. Б. Метод продолжения решения по параметру и наилучшая параметризация. М.: Эдиториал УРСС, 1999.
[7] Фадеев С. И., Покровская С. А., Березин А. Ю., Гайнова И. А. Пакет программ STEP для численного исследования нелинейных уравнений и автономных систем общего вида. Новосибирск: НГУ, 1998.
[8] Корнев В. М., Коробейников С. Н. О реализуемости симметрии решения в задаче осесимметричного деформирования цилиндрической оболочки при продольном
сжатии / / Аналитические и численные методы решения краевых задач пластичности и вязкоупругости. Свердловск: УНЦ АН СССР, 1986. С. 70-77.
[9] Аннин Б. Д., Алехин В. В., КОРОБЕЙНИКОВ С.Н. Определение предельных состояний упругопластических тел // ПМТФ. 2000. Т. 41, №5. С. 196-204.
[10] БУРАГО Н.Г., Кукуджлнов В. Н. Численный метод решения геометрически нелинейных осесимметричных задач для упругопластических оболочек вращения // Строительная механика и расчет сооружений. 1976. №5. С. 44-49.
[11] Kornev V. M., Korobeinikov S. N. On solution symmetry for axisymmetric deformation of axially compressed cylindrical shell // Innovative Numer. Methods in Eng.: Proc. 4th Intern. Symp. / Eds R.P. Shaw et al. Berlin: Springer-Verlag, 1986. P. 579-584.
[12] Коробейников С. Н. Нелинейное деформирование твердых тел. Новосибирск: Изд-во СО РАН, 2000.
[13] STEPHENS W. B. Computer program for static and dynamic axisymmetric nonlinear response of symmetrically loaded orthotropic shell of revolution. 1970 (Report / NASA; TN D-6158).
[14] Hill R. On uniqueness and stability in the theory of finite elastic strain //J. Mech. Phys. Solids. 1957. Vol. 5, No. 4. P. 229-241; О единственности и устойчивости в теории конечных упругих деформаций // Механика (сб. переводов). 1958. №3(49). С. 53-65.
[15] Коробейников С. Н. Достаточное условие потери устойчивости квазистатических движений деформируемых тел // Симметрия и дифференциальные уравнения: Тр. Междунар. конф. / СО РАН. Ин-т вычислительного моделирования. Красноярск, 2000. С. 125-127.
[16] Хилл Р. Бифуркация и единственность в нелинейной механике сплошной среды // Проблемы механики сплошной среды: Сб. науч. тр. М.: АН СССР, 1961. С. 448-457.
[17] Hill R. A general theory of uniqueness and stability in elastic-plastic solids // J. Mech. Phys. Solids. 1958. Vol. 6, No. 3. P. 236-249; Общая теория единственности и устойчивости для упругопластических тел // Механика (сб. переводов). 1958. №6(52). С. 81-95.
[18] Клюшников В. Д. Устойчивость упругопластических систем. М.: Наука, 1980.
[19] STRANG G. Linear Algebra and Its Applications. N. Y.: Acad. Press, 1976; Линейная алгебра и ее применения. М.: Мир, 1980.
[20] BuSHNELL D. Effect of cold bending and welding on buckling of ring-stiffened cylinders // Computers & Structures. 1980. Vol. 12. P. 291-307.
[21] ФАДДЕЕВ Д. К., ФАДДЕЕВА В. Н. Вычислительные методы линейной алгебры. М.; Л.: Физматгиз, 1963.
[22] БАЕВ Л. В., Коробейников С. Н. Выпучивание круговой цилиндрической оболочки за пределом упругости при действии осевой силы и бокового давления // Динамика сплошной среды / АН СССР. Сиб. отд-ние. Ин-т гидродинамики. Новосибирск, 1982. Вып. 55. С. 9-31.
[23] Tvergaard v., Needleman A. On the development of localized buckling patterns // Collapse: the Buckling of Structures in Theory and Practice: Proc. Int. Symp. / Eds: J. M. T Thompson, G.W. Hunt. Cambridge: Cambridge Univ. Press, 1983. P. 11-24; О явлении локализованного выпучивания // Потеря устойчивости и выпучивание конструкций: теория и практика. М.: Наука. 1991. С. 11-24.
[24] Hutchinson J.W. Plastic buckling // Advances in Applied Mechanics. Vol. 14 / C. S. Yih (Ed.) - N. Y.: Acad. Press, 1974. P. 67-144.
[25] Sewell M. J. The static perturbation technique in buckling problems //J. Mech. Phys. Solids. 1965. Vol. 13, No. 4. P. 247-265.
[26] Hill R. Eigenmodal deformations in elastic/plastic continua //J. Mech. Phys. Solids. 1967. Vol. 15, No. 6. P. 371-386.
[27] Коробейников С. Н. Применение метода конечных элементов к решению нелинейных задач по деформированию и потере устойчивости атомных решеток. Новосибирск, 1997 (Препр. / РАН. Сиб. отд-ние. Ин-т гидродинамики, №1-97).
[28] Кирсанов М. Н. Модельное представление потери устойчивости при неоднородном докритическом состоянии. М., 1982 (деп. ВИНИТИ. №303).
[29] Вольмир А. С. Устойчивость деформируемых систем. М.: Наука, 1967.
[30] Григолюк Э. И., Кабанов В. В. Устойчивость оболочек. М.: Наука, 1978.
[31] Geckeler J.W. Plastisches knicken der wandung von hohlzylindern und einige andere faltungserscheinungen an schalen und blechen // ZAMM. 1928. Vol. 8, No. 5. P. 341-352.
[32] Johnson W., Soden P. D., Al-Hassani S. T. S. Inextensional collapse of thin-walled tubes under axial compression //J. Strain Analysis. 1977. Vol. 12, No. 4. P. 317-330.
[33] Зубчанинов В. Г. Экспериментальное исследование процессов потери устойчивости цилиндрических оболочек при осевом сжатии // Инженерный журнал. 1965. Т. 5, №3. С. 583-586.
[34] Тетерс Г. А. Сложное нагружение и устойчивость оболочек из полимерных материалов. Рига: Зинатне, 1969.
Поступила в редакцию 26 апреля 2000 г.
