ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ЗАКРИТИЧЕСКОГО НЕЛИНЕЙНОГО ДЕФОРМИРОВАНИЯ ОСЕСИММЕТРИЧНЫХ МЕМБРАН Текст научной статьи по специальности «Математика»

УДК:539.38: 621.01: 004.7

БОТ: 10.18698/2309-3684-2020-1-6487

Численное моделирование закритического нелинейного деформирования осесимметричных мембран

© С. А. Подкопаев МГТУ им. Н.Э. Баумана, Москва, 105005, Россия

Рассмотрены теоретические основы нелинейного деформирования тонких осесимметричных оболочек. Представлены эксплуатационные характеристики мембран в различных коммутационных устройствах, клапанах и датчиков давления. Рассмотрены типы нелинейного поведения закритического поведения мембран. Представлена математическая модель для описания процесса нелинейного деформирования осесимметричных оболочек, метод дискретного продолжения по параметру и прием «смены подпространства управляющих параметров». На примере шарнирно-опертой сферической оболочки выполнено исследование закритического поведения. Выбрана рациональная математическая модель для описания нелинейного деформирования хлопающих симметричных оболочек. Разработан и реализован в виде авторской программы численный алгоритм исследования процессов нелинейного деформирования многопараметрических систем.

Ключевые слова: нелинейное деформирование, тонкостенная осесимметричная оболочка, мембрана, закритическое поведение, дискретное переключение, продолжение по параметру, смена подпространства параметров

Введение. Разнообразные элементы конструкций и технических устройств, выполненные в формы осесимметричных оболочек (мембран), широко применяются и используют в различных отраслях промышленности. В настоящее время, в связи цифровой промышленной революцией (Индустрия 4.0), особенно массовое распространение получили различные коммутационные устройства, предохранители и клапаны, применяемые в промышленном интернете вещей [2, 6, 11, 12]. Под мембраной (рис. 1а) будем принимать тонкостенную осесимметричную оболочку, которая, под действием внешней нагрузки может скачкообразно менять свой прогиб. Данный способ потери устойчивости, не сопровождающийся разрушением оболочки, будем называть прощелкиванием [14, 15, 16].

Эксплуатационные характеристики мембран, изготовленные в виде оболочек. Мембраны в форме оболочке находит применение в различных коммутационных устройствах (рис. 2а), клапанах, сигнализаторах и датчиках давления (рис. 2б). Осесимметричная оболочка в виде мембраны также применяется в конструкции топливных баков авиастроении и кораблестроении в качестве вытеснительной мембраны (рис. 2в).

Главным эксплуатационным элементом мембраны является упругая характеристика, отражающая зависимость между перемеще-

нием контрольной точки мембраны и изменением внешней нагрузки (рис. 16) [13, 17].

Показанные на рис. 16 участки графика АВ, СО устойчивые участки упругой характеристики мембраны. Под воздействием внешней нагрузки при достижении первой экстремальной точки (экстремума) (точка В, рис. 1б), соответствующей величине критического давления р , мембрана, минуя неустойчивый участок графика ВС, скачкообразно меняет прогиб и процесс деформирования продолжиться по участку СО упругой характеристики. При уменьшении (снятии) внешней нагрузки с оболочки будет \'0 происходить обратное скачкообразное изменение прогиба, соответствующее второму критическому

В

Г

1 Ркрх '■. / •-.У

1А Ркр2

Рис. 1. Мембрана

а — пример мембраны; б — вид упругой характеристики мембраны

давлению ркр^. Следует отметить,

что после первоначальной нагрузки (первый скачок) напряжения в мембране будут оставаться упругими.

Рис. 2. Применение мембран в форме оболочек

а — коммутационное устройство (кнопка); б — клапан (1 — мембрана; 2 — корпус); в —топливный бак (В — топливо; Г — газ; М — вытесняющая мембрана)

а

в

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

отличных от р и р. Давление, при котором оболочка будет прощелкивать, будем обозначать как рхл. Значение давления рхл лежит в интервале между верхним рк^ и нижним р критическими давлениями.

Типы исследования закритического поведения мембран.

Можно выделить два основных подхода к исследованию закритиче-ского поведения нелинейных механических систем.

Первый подход - классический. Он заключается в нахождении критических значений величины внешней нагрузки, при которых существуют смежные формы равновесия конструкции [1, 5, 9].

Второй подход заключается в непосредственном построении самой поверхности равновесных состояний в пространстве параметров системы. Достоинством данного подхода является возможность детального исследования закритического поведения оболочки, которое зачастую и характеризует эксплуатационные свойства упругого элемента [3, 4, 6, 7, 23, 24, 25].

В статье авторами предлагается и используется многопараметрический подход к исследованию закритического нелинейного процесса деформирования прохлопывающих оболочек (мембран).

Предлагаемый подход использует идею, при которой происходит процесс «погружения» конкретной задачи в многопараметрическое семейство подобных задач.

Такой подход является удобным инструментом для проектирования и конструирования технических конструкций, позволяющий получить не одно частное решение, а решение сразу целого семейства задач прохлопывающих мембран, в зависимости от внешних параметров системы.

Для построения поверхности равновесных состояний мембраны был использован дискретный метод продолжения по параметру одновременно с приемом «смены подпространства управляющих параметров». Иначе говоря, сложную многопараметрическую задачу можно рассматривать, как множество однопараметрических задач, каждая из которых имеет свой дискретно изменяющийся параметр. Для перехода между однопараметрическими задачами был использован прием смены подпространства управляющих параметров [6, 7, 10, 18, 19].

Задача исследования закритического поведения оболочки фактически сводится к анализу поверхности равновесных состояний. Для некоторых значений параметров системы может наблюдаться многозначность решений, наличие особых точек и изолированных решений. Отличительная особенность данного подхода заключается в том, что вместо решения бифуркационной задачи, предлагается численно

исследовать графики перестроек решения при переходе через особую точку [4, 6, 7, 21, 22].

Для построения графика перестроек и обхода окрестностей бифуркационных точек применяется прием «смены подпространства внешних параметров».

Использование метода конечных элементов (МКЭ) при исследовании закритического поведения оболочек связано с рядом трудностей, таких как с отсутствием возможности движения по параметрам геометрии, по толщине или по радиусу кривизны. Иначе говоря, варьирование любыми геометрическими параметрами конструкции повлечет за собой многократное перестроение конечно-элементной сетки, что в свою очередь связано со значительными требуемыми вычислительными мощностями.

В частности, МКЭ не позволяет отыскать изолированные решения на кривой равновесных состояний, а также смоделировать процесс зарождения на упругой характеристике так называемых «петель» [18, 20].

Математическая модель для описания процесса нелинейного деформирования осесимметричных оболочек. В статье были использованы соотношения теории тонких осесимметричных оболочек, модернизированные для последующего применения численного алгоритма. При описании геометрии мембраны (оболочки вращения) в качестве отсчетной поверхности была использована срединная поверхность оболочки (рис. 3а) Точка В0 на срединной поверхности характеризуется двумя координатами: углом р, определяющим положение меридионального сечения, и дуговой координатой £0, направленной вдоль меридиана. При осесимметричной деформации все меридиональные сечения оболочки равноправны и, следовательно, существенной является только координата £0.

Геометрия меридионального сечения мембраны (осесимметрич-ной оболочки) (рис. 3б) в исходном недеформированном состоянии была задана в параметрическом виде:

Х0 = Х0(Б0), ¥0 = 70(^0), (1)

где £0 — независимая координата, которая отсчитывается от предварительно зафиксированной точки меридиана А до текущей точки В меридиана; Х,У — декартовая система координат, в которой ось X направлена по радиусу оболочки, а ось У совпадает с осью вращения оболочки; Х0, у — декартовы координаты текущей точки меридиана до деформирования.

Обозначим текущий угол наклона касательной к меридиану в не-деформированном состоянии как щ .

Запишем геометрические соотношения:

Рис. 3. К выводу основных соотношений осесимметричной оболочки а — срединная поверхность оболочки; б — геометрия меридионального сечения мембраны (осесимметричной оболочки)

Параметры, относящиеся к недеформированному состоянию, будут обозначаться с использованием нижнего индекса "0".

После деформации точка В0, принадлежащая срединной поверхности оболочки, переместится в новое пространственное положение В, характеризующееся координатами X, У и Б, соответственно (рис. 3б). Тогда геометрические соотношения для деформированного состояния запишутся в виде:

с1Х сСУ .

-= сощ, — = Бтщ. (3)

СБ СБ

Перемещения по оси X - и и перемещения по оси У - V точки В и угол поворота нормали в определяются следующими соотношениями:

и = X-Х0, V = У-У0, в = щ-щ . (4)

Для обозначения величин, соответствующих меридиональному и окружному направлениям, используется нижний индекс «1» и «2» соответственно.

Выражения для главных радиусов кривизны в исходном и деформированном состояниях оболочки имеют вид:

1 _ dy0 _ 1 _ dy 1 _siny0 1 _siny Ло dS0' P1 dS' P20 X0 ' P2 X

Полные изменения главных кривизн обозначим следующим образом:

йщ, _ siny sin^0 (еЛ

Kj — , Л2 — . (О/

dS dS0 X Xo

Полное изменение кривизны в меридиональном направлении складывается из изменения кривизн вследствие поворота нормали (изгиба) и удлинения меридиана:

Ccompl Сbend ^ Сgeom .

Тогда:

= _ J dyL W dyL_ dyL \ ёщ^

С bend С compl С geom ^^ ^^ то 1Г, TO TO . V7/

^ dS dSo j

dS dS

Введем обозначение:

. (8)

10 dS dS

Аналогичным образом изменение кривизны в окружном направлении:

20 х х ■ ( )

Линейная деформация элемента срединной поверхности в меридиональном направлении для текущего состояния имеет вид:

,10=dS_dSо. (Ш)

10 dS0

Линейная деформация элемента срединной поверхности в окружном направлении для текущего состояния имеет вид:

_(х0 + и)'У_х0 'У _ и

£20 ~ V т ~ V ' (11)

Х0 'У Х0

Преобразуем выражение (8) к следующему виду:

1

«10 =

г сСщ Сщ0 л

(1 + ^10)

СБа СБ(

(12)

о у

Используя соотношения (1-4), (10), (12), получим следующие геометрические соотношения уравнений осесимметричной тонкостенной оболочки:

йи

СБ0 СV

СБ0 Сщ

= (1 + 810) • сощ - со8щ0;

= (1 + 810) • $тщ- 8тщ0;

(13)

СБа

= (1 + 810) «10 +

Що СБа

Напряжено-деформированное состояние осесимметричной оболочки. Эквидистантная поверхность - поверхность, равноудаленная от срединной поверхности на расстоянии . Расстояние отсчитывается по направлению внешней нормали к недеформиро-ванной срединной поверхности. Причем нормальная координата текущей материальной точки изменяется в следующих пределах:

|(Бо) |(Бо),

(14)

где к(Б0) — толщина оболочки, в общем случае зависящая от дуговой координаты Б .

С учетом предположения о тонкостенности оболочки:

тах[/г(Б0)]«тт(|/?10|,|/?20|). (15)

Примем кинематическую гипотезу Киргофа-Лява, в соответствии с которой материальный отрезок нормальный к срединной поверхности оболочки до деформации и после деформации остается нормальным к срединной поверхности деформированной оболочки. При этом длина этого отрезка не изменяется .

Тогда для линейных деформаций произвольного элемента в малой окрестности точки с координатами Б0 и справедливы следующие соотношения:

81| (Б0, ¿0 ) = 810 (Б0 ) + ¿0 • «10 (Бо ), 82| (Б0, ¿0 ) = 820 (Б0 ) + ¿0 • «20 (Б0 ).

(16)

Следствием выражений (16) является линейное распределение деформаций по толщине оболочки, и то, что деформация эквидистан-той поверхности будет полностью определена, если будут известны мембранные и изгибные деформации срединной поверхности.

Примем статическую гипотезу Киргофа-Лява, согласно которой можно пренебречь нормальными напряжениями <г3, действующими

на площадках параллельных срединной поверхности оболочки.

В соответствии с принятой гипотезой, оболочка испытывает плоское напряженное состояние. Запишем закон Гука для плоского напряженного состояния:

Е

<1 (17)

Внутренние силы, возникающие в сечениях оболочки, приведем к ее срединной поверхности (рис. 4):

N1 =/\ М =/\ <1

"I "I (18)

N2 = \ \ <2^ М2 =\ \ <2 • ^.

2 2 При интегрировании в формулах (18), в силу тонкостенности

оболочки, пренебрегаются величины порядка < по сравнению с

Р

единицей. Последовательно подставляя в выражение (18) соотношения (16-17), приходим к следующим выражениям для силовых факторов:

^ = С1 • £10 + С2 • ^20; М1 = С3 • «10 + С4 • «20; N2 = С1 £20 + С2 £10; М2 = С3 • «20 + С4 «10.

(19)

Рис. 4. Силовые факторы, действующие в срединной поверхности элемента оболочки

Постоянные коэффициенты определяются следующим образом:

_ Е п _ /лЕк С1 — , 2 ' С2 — - ;

С3 —

1 -л

Ек3

1 -л

; С4 — ■

лЕк

(20)

12 (1 -л2 У 4 12(1 -л2 )• Введем следующие обозначения:

5 —

Ек 1 -л2

Б —

Ек3

12 (1 -л2 )'

(21)

где 5 — мембранная жесткость оболочки; Б — цилиндрическая жесткость оболочки.

Вывод уравнений равновесия сил и моментов. Меридиональная N и поперечная 0 — силы, связанные с горизонтальной и и вертикальной V составляющими внутреннего усилия при помощи матрицы поворота, зависящей от угла щ (рис. 5).

[N1'

101,

СОБЩ БШЩ

- ътщ собщ

и 1

V I

(22)

Рис. 5. Силы, действующие на элемент срединной поверхности оболочки (к выводу уравнений равновесия осесимметричной оболочки)

Запишем уравнение равновесия сил элемента оболочки. Уравнения равновесия запишем в глобальной системе координат X, У. Для этого спроектируем все силы, действующие на элемент срединной поверхности оболочки, на направления X и У (рис. 5).

й(X■V■ йф) + ч ■ X■ йу-йБ - 0;

йф (23)

й(X■и■ йф)-2ЩйБ■ -ф + ди ■ X■ йф-йБ - 0.

После преобразования (23) получаем следующие уравнения:

й (X ■V)

йБ0 й (X и) йБп

--(1■X;

--(1 + X

N2..

^ Чи VЛ у

(24)

Запишем уравнение равновесия моментов всех внешних и внутренних сил, приложенных к элементу оболочки, относительно оси I. Для этого векторы моментов необходимо спроектировать на направление I (рис. 6б).

йй

й(М1 ■ X) ■ йф-2М2 ■ йБ ■ — + ( ■ X■ йф■ йБ - 0. (25)

б

Рис. 6. Моменты, приложенные к элементу оболочки, относительно оси I (к выводу уравнений равновесия момента

Рассматривая рис. 6а, выразим приращение угла йй следующим образом:

Й? й- й ф. (26)

Учитывая соотношения (10) и (22), преобразуем уравнение (25) к следующему виду:

й (М ■ X)

- (1 + %) ■ X

М2---ъи ■ V ■ соъщ

X

- 0. (27)

а

Принимая во внимание соотношения (3-4), преобразуем уравнения (24) и (27) к следующему виду:

Л (V )

Л (и )

ЛБ0

= "(1 + ^10 )

совщ

Х0 + и

V + я,

= "(! + ею )

сов щ гг N

-и "■ 2

Х0 + и Х0 + и

Яи

Л (м)

= "(1 + £10 )

(М " М 2 )■

сощ Х +и

" и ■ в1п щ + V ■ сов щ

(28)

где яи и — интенсивности распределенной нагрузки.

ЯУ = Яп сощ + Яг втщ; Яи ="Яп 81ПЩ + Яг совщ.

(29)

Двухточечная краевая задача для системы нелинейных дифференциальных уравнений шестого порядка. Соотношения (13), (28) образуют основную систему нелинейных дифференциальных уравнений.

Неизвестные, производные которых входят в уравнения, будем называть основными неизвестными. Тогда вектор основных неизвестных в текущем сечения оболочки запишется в следующем виде:

{ X }Т ={и, у,щ,и, V, М}.

(30)

Все остальные неизвестные будем называть вспомогательными. Для решения основной системы вспомогательные неизвестные необходимо выразить через основные неизвестные.

Используя соотношения (19), (20), получаем следующие формулы для нахождения вспомогательных величин:

1

е10 = —(и совщ + V в1п щ)" ^—;

В

_М1 «10 = " ^

Х0 + и V Х0 У

Х

81И Щ 81И Щ0

Л

V Х0 + и

Х

0 У

и

N = м(и совщ + V в1пщ) + ЕИ-;

М2 = /иМ1

ЕИ3

12

Х0 + и

Х

Х

81П Щ 81П Щ0

Л

V Х0 + и

Х

0 У

(31)

<

Дополняя систему (13), (28) соответствующими краевыми условиями сводим проблему исследования потери устойчивости мембраны к двухточечной краевой задаче для системы нелинейных дифференциальных уравнений шестого порядка в обыкновенных производных, которую представим в виде:

йX

йБг,

Р (Б,, Z, X, ();

С0 (X, ()- 0; О, (X, ()- 0,

(32)

где ( — векторный параметр внешней нагрузки; Z — вектор вспомогательных неизвестных, не входящих под знак производной.

Переход от многопараметрической задачи к множеству одно-параметрических задач. Метод дискретного продолжения по параметру. Для исследования сложных многопараметрических процессов деформирования семейство нелинейных краевых задач в обыкновенных производных было сведено к многопараметрическому семейству систем нелинейных уравнений вида:

р ({X},{X2})- 0. (33)

Система (33) порядка т зависит от вектора «внутренних» параметров {размерностью т х1 и вектора «внешних» параметров

{X2} размерностью п х1. В состав вектора {X1} входят параметры,

характеризующие состояние рассматриваемой системы, прежде всего, это обобщенные перемещения и обобщенные внутренние силовые факторы. Внешние параметры или по-другому параметры «управления», входящие в состав вектора {X2} являются переменными, и, как

правило, связаны с геометрическими размерами конструкции, свойствами материала, условиями закрепления, внешней нагрузкой и т.п.

Совокупность всех решений (33) можно рассматривать как некоторую поверхность (гиперповерхность) равновесных состояний, построенную в евклидовом пространстве размерности Ят+п. Построение такой гиперповерхности является очень трудоемким процессом, поэтому, зачастую решение сложной многопараметрической задачи сводят к решению семейства однопараметрических задач, каждая из которых имеет свой, дискретно изменяющийся параметр. Такой подход позволяет построить сечения поверхности равновесных состояний, что значительно упрощает исследуемую задачу и дает возмож-

ность проанализировать влияние «внешних» параметров на характер деформирования рассматриваемого элемента. В этом случае вектор «внешних» параметров {Х2} определяется через один независимый скалярный параметр я •

{ Х2 } = Я. (34)

Для исследования задачи потери устойчивости хлопающей мембраны, нагруженной внешним давлением, применялся метод продолжения по параметру, основная идея которого, заключается в том, что любую сложную многопараметрическую задачу можно рассматривать, как множество однопараметрических задач, каждая из которых имеет свой, дискретно изменяющийся параметр при фиксированных значениях всех остальных внешних параметрах.

В качестве варьируемого параметра Я может быть взята внешняя нагрузка, геометрическая характеристика системы, свойство материала, условия закрепления и т.п.

Для дальнейшей алгоритмизации удобно считать параметр Я равноправным со всеми внутренними параметрами системы, то есть с компонентами вектора состояния {Х} системы (32).

Введем расширенный вектор основных неизвестных следующим образом:

{Х _ ехг}={{ Х}, я}. (35)

Тогда система (33) для однопараметрического процесса запишется следующим образом:

г ({Х _ ехгг}) = 0. (36)

Двухточечная процедура по схеме «предиктор-корректор».

При реализации процедуры дискретного продолжения по параметру использовалась двухэтапная процедура по схеме «предиктор-корректор». На этапе «предиктор» на основании предыстории при помощи экстраполяции решения осуществляется предсказание

начального вектора {Х^г} = {{Хк }, як } для нового значения параметра як ■

Верхний индекс к соответствует номеру итерации по параметру. Для возможности проведения экстраполяции необходимо сохранять информацию о ранее полученных решениях.

Однако исключением является случай начала итерационного процесса. В такой ситуации необходим так называемый этап «разгонки». На первом шаге по параметру экстраполяция начального

приближения не проводится, а в качестве начального приближения, используется некоторое «опорное» решение. Применительно к рассматриваемой задаче деформирования осесимметричной оболочки вращения, в качестве «опорного» решения используется тривиальное решение, соответствующее нулевой внешней нагрузке. Также в качестве «опорного» решения можно использовать решение задачи в линейной постановке, либо воспользоваться решением метода конечных элементов.

На втором шаге применяется линейная экстраполяция по двум точкам. И только потом, получив решение в трех точках равновесной кривой, процедура выводится на стандартный режим квадратичной экстраполяции.

На этапе «корректор» происходит уточнение начального приближения {X0} с помощью итерационного метода. В качестве итерационного метода используется метод Ньютона и его модификации.

Модифицированный метод Ньютона-Рафсона. На этапе «корректор» итерационным методом Ньютона решается нелинейная двухточечная краевая задача, посредством многократного решения задачи Коши, и итерационного уточнения начального вектора для системы (36).

Получив на этапе «предиктор» начальное приближение вектора

основных неизвестных

{X0* 1

4=0

соответствующего текущей величине

параметра , подставляем его в систему нелинейных дифференциальных уравнений (32). Решив задачу Коши с приближенным вектором основных неизвестных, получаем некоторую невязку решения (36). Так как все дальнейшие выкладки в данном разделе соответствуют одному шагу по параметру к , обозначим вектор основных

неизвестных {X* } как {Xп }, где верхний индекс будет отвечать

' 4=0 ^ '

номеру итерации в методе Ньютона.

При п = 0 вектор {Xп} имеет следующий вид:

{X0}

(37)

где х, X, х3 — известные компоненты вектора | X01; х° —

компоненты, найденные на этапе «предиктор», требующие уточне ния.

После подстановки (37) в систему (32) получим: й { X0 }

йз

Е (з,{ X0}, а ) = |г0

(38)

где {г0} — вектор невязки нулевого приближения или номинальная невязка.

На рис. 7 приведена геометрическая интерпретация модифицированного метода Ньютона-Рафсона (метода одной касательной).

В точке начала интегрирования задается начальное приближение вектора основных

{г}

{г °) {г0}' {г'}

неизвестных

{X0}, которому

{X}

И И И

Рис. 7. Геометрическая интерпретация метода Ньютона (метода одной касательной)

соответствует вектор невязки {г0 } . Затем поочередно каждой неизвестной компонент

X (г = 4,5,6) вектора {X0} задается малое приращение (возмущение) Лхг, тем самым получаем три новых «возмущен-

ных» вектора решения {X0} . Проведем через точку С касательную к графику исследуемой функции. Точка пересечения касательной с осью абсцисс берется в качестве следующего приближения {X'}.

Для нового приближения {X'} все действия аналогичные: определяется новая невязка решения {г'}, ей соответствует новая точка на

кривой, в этой точке проводится касательная, точка пересечения касательной с осью абсцисс определяет следующее приближение. Итерационный процесс считается сошедшимся, если конечная невязка будет меньше некоторого заранее выбранного значения 8. Определим тангенс угла наклона касательной:

Дг,

1ап а = ■

Дх

(39)

где Дт — приращение невязки \ -ого неизвестного компонента вектора {X0|; Дху. — приращение ] -ого неизвестного компонента

вектора {X0|.

Получим выражение для нового приближения | X11 (рис. 7):

i х1Н х {'"Í

tana =4—^-^{ДХ } = {ДХ} 1 j

tana

(40)

(41)

Подставив (41) в выражение (40) и обобщив для многомерного случая, получим основное рекуррентное соотношение метода Нью-тона-Рафсона:

{Хп+1 } = { Хп}-[ J0 ]_1 {rn},

(42)

где |Х"+1| — вектор решения (п +1) -ой итерации; |Хп | — вектор решения п -ой итерации; |гп | — вектор невязки п -ой итерации;

30 — матрица Якоби, вычисленная в нулевом приближении.

При численном счете аналогом матрицы Якоби является матрица Гато, которая имеет следующий вид:

J0

Дг Дг Дг

Дх1 Дх2 Дхз

ДГ2 _ДГ2_ ДГ2

Дх1 Дх2 Дх2

Дгз Дгз Дгз

Дх1 Дх2 Дхз

(43)

Компоненты матрицы Гато вычисляются следующим образом:

Дт = т ({X0О)-т ({X0}). (44)

Отличие модифицированного метода Ньютона-Рафсона (метод одной касательной) от классического состоит в том, чтобы вычислять матрицу Якоби лишь один раз, в точке начального приближения, а затем использовать это значение на каждой последующей итерации.

Прием «смены подпространства управляющих параметров».

Для обхода точек бифуркации или точек ветвления, возникающих на кривой равновесных состояний, используется прием численного счета, получивший название «прием смены пространства управляющих параметров» [4, 6, 7, 11, 18].

Основная идея приема заключается в том, что при подходе к окрестности предполагаемой точки бифуркации следует перейти к новой системе, для которой на кривой равновесных состояний уже не будет никаких бифуркационных точек или точек ветвления (наличие предельных точек допускается). Фактически переход к новой системе можно интерпретировать как решение задачи с немного измененной конфигурацией, которая принадлежит целому семейству подобных задач. После прохождения критического участка (окрестности бифуркационной точки) можно совершить обратный переход и продолжить решение задачи с исходной конфигурацией.

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

Тогда однопараметрическая система должна быть записана следующим образом:

г ({X!}, д1г д2 ) = г ({Хех(г }) = 0, (45)

где {Хех(г} = {{Х1},д2} — расширенный вектор; {Х1} — вектор основных неизвестных, имеет размерность т х1; ^, д2 — два параметра управления.

Система (44) имеет размерность т + 2.

Исследование закритического поведения шарнирно-опертой сферической оболочки. Рассмотрим многопараметрический подход и прием «смены подпространства управляющих параметров» на примере исследования закритического поведения сферической оболочки, шарнирно опертой по внешнему контуру и нагруженной внешним давлением (рис. 8). В таблице 1 приведены исходные данные.

1/

Рис. 8 Схема сферической оболочки, шарнирно опертой по внешнему контуру

Таблица 1

Исходные данные

Радиус опорной поверхности а = 2,8 мм

Толщина оболочки Н = 0,05 мм

Исходный радиус кривизны меридиана Я] = 32 мм

Модуль упругости Е = 13104 МПа

Коэффициент Пуассона ц = 0,3

Способ закрепления края оболочки Шарнир

Многопараметрический подход использует стратегию последовательного исследования однопараметрических семейств размерности Ят+1, каждая из которых принадлежит многопараметрическому

семейству задач размерности Ят+п. Каждой однопараметрической задаче соответствует свой дискретно изменяющийся внешний параметр, при этом остальные (п -1) управляющие (внешние) параметры

имеют фиксированные значения.

В примере исследование закритического поведения сферической оболочки проводилось в пространстве: прогиб оболочки в центре (V), радиус кривизны (Я), внешнее давление (р). Радиус кривизны

и внешнее давление являются управляющими параметрами, а прогиб оболочки в центре — одна из компонент вектора внутренних параметров, характеризующих текущее состояние конструкции.

На рис. 9 представлена проекция поверхности равновесных состояний в пространство: V, Я, р. Данная поверхность построена с использованием метода продолжения по параметру в связке с приемом «смены подпространства управляющих параметров». На рисунке видно, что с увеличением радиуса кривизны проекция поверхность равновесных состояний усложняется, и при некотором значении параметра кривизны появляется отдельная поверхность изолированных решений, которая впоследствии сливается с основной поверхностью.

Продемонстрируем прием «смены подпространства управляющих параметров», для этого рассмотрим сечения проекции поверхности равновесных состояний для значений радиуса кривизны Я , лежащего в пределах от 32 мм до 36 мм (рис. 10). Траектории 1 соответствует радиусу кривизны равному 35,5 мм, а траектория 3 соответствует радиусу кривизны равному 32 мм. Наблюдая за качественным изменением траекторий при монотонном изменении параметра кривизны, можно предположить, что существует некоторая особая точка — точка бифуркации, соответствующая некоторому критическому значению параметра кривизны Я

(33,5185 < Якр < 33,5188). На рис. 10 сечение, соответствующее критическому значению параметра кривизны, изображено заштрихованной плоскостью.

р, МПа

0,35 0,3 0,25 0,2 0,15 0,1 0,05 0

-0,05

-0,1 44

40 38 Я, мм

36 0

0,15 V, мм

0,2 0,25°

0,35

Рис. 9. Проекция поверхности равновесных состояний в пространстве: внешнее давление (р); радиус кривизны (Я); прогиб в центре (V)

р, МПа

0,25

0,2

0,15

0,1

0,05

0

0,05

-0,1

0

0,05

35

35,5

36

Я, мм

Рис. 10 Сечения проекции поверхности равновесных состояний в пространстве: внешнее давление (р), радиус кривизны (Я), прогиб в центре (V)

Приближаясь к критическому значению параметра кривизны Я снизу, траектории, лежащие в плоскости прогиб давление, будут

иметь вид аналогичный траектории 3, а если приближаться сверху, то аналогичный траектории 1.

Получить точное значение критического параметра оказывается практически невозможным. По мере приближения к окрестности бифуркационной точки наблюдается ухудшение сходимости численного решения: может наблюдаться самопроизвольный переход на другую ветвь решения, либо разворот назад.

В данной статье в рамках многопараметрического подхода предлагается альтернативная стратегия, позволяющая избежать непосредственного решения задачи ветвления, но при этом получить всю необходимую информации о поведении системы в окрестности бифуркационной точки. Иначе говоря, можно обойти особую точку с разных сторон и построить так называемую картину перестройки.

С помощью приема смены подпространства управляющих параметров можно в нужный момент ответвится от траектории 3, соответствующей однопараметрической задаче с переменным давлением (р = var) и фиксированной кривизной (Я = 32 мм), и начать движение вдоль траектории 4 (рис. 10), которая соответствует другой однопараметрической задаче уже с переменной кривизной (Я = var), но

с фиксированным внешним давлением (р = 0,084 МПа).

Фактически, прием смены подпространства управляющих параметров является переключателем между разными однопараметриче-скими задачами.

Достигнув заданного значения параметра радиуса кривизны, например, Я = 35,5 мм, можно повторно ответвиться от траектории 4, и таким образом попасть на траекторию 5 (рис. 10), соответствующую изолированному решению и не имеющую никаких общих точек с основной траекторией 1.

Продолжив движение вдоль траектории 4, можно определить «глубину» проникновения изолированного решения. Как видно из полученных результатов, траектория 4 также имеет предельную точку (в данном случае Япр = 35,8 мм). Следовательно, при значениях

параметра радиуса кривизны больших, чем Япр изолированное решение, аналогичное траекториям 5 и 6, отсутствует. Для преодоления предельной точки на траектории 4 используется прием смены параметра продолжения.

Продолжая движение по траектории 4, минуя предельную точку, можно снова вернуться на исходную траекторию 3. Таким образом, осуществляется обход окрестности бифуркационной точки.

Данный алгоритм был реализован в авторской программе и для подтверждения достоверности результатов численного решения рассматриваемая задача была решена с использованием метода конечных элементов в программном комплексе ANSYS. Результаты автор-

ской программы с высокой точностью совпали с результатами программного комплекса ANSYS, что свидетельствует о достоверности представленных результатов.

Выводы. В статье выполнен обзор существующих подходов и методик исследования закритического поведения осесимметричных сферических оболочек.

1. Выбрана рациональная математическая модель для описания нелинейного деформирования хлопающих симметричных оболочек.

2. Разработан и реализован в виде авторской программы численный алгоритм исследования процессов нелинейного деформирования многопараметрических систем. В рамках алгоритма была разработана методика расчета однопараметрических задач с использованием дискретного метода продолжения по параметру, был реализован прием «смены параметра продолжения» для преодоления особых предельных точек кривой равновесных состояний.

ЛИТЕРАТУРА

[1] Алфутов Н.А. Основы расчета на устойчивость упругих систем. Библиотека расчетчика. Москва, Машиностроение, 1978, 312 с.

[2] Андреева Л.Е. Упругие элементы приборов: учебное пособие. Москва, Машиностроение, 1982, 456 с.

[3] Бидерман В.Л. Механика тонкостенных конструкций. Статика. Библиотека расчетчика. Москва, Машиностроение, 1977, 488 с.

[4] Валишвили Н.В. Методы расчета оболочек вращения на ЭЦВМ. Москва, Машиностроение, 1976, 278 с.

[5] Вольмир А.С. Устойчивость деформируемых систем. Москва, Физматгиз, 1967, 984 с.

[6] Гаврюшин С.С. Разработка методов расчета и проектирования упругих оболочечных конструкций приборных устройств: диссертация. Москва, 1994, 316 с.

[7] Гаврюшин С.С., Барышникова О.О., Борискин О.Ф. Численные методы в динамике и прочности машин. Москва, Издательство МГТУ им. Н.Э. Баумана, 2012, 492 с.

[8] Григолюк Э.И., Лопаницын Е.А. Конечные прогибы, устойчивость и за-критическое поведение тонких пологих оболочек. Москва, МГТУ «МА-МИ», 2004, 162 с.

[9] Григолюк Э.И., Кабанов В.В. Устойчивость оболочек. Москва, Наука, 1978, 360 с.

[10] Ильгамов М.А. О научном наследии Х.М. Муштари. Труды Математического центра им. Н.И. Лобачевского, 2010, т .42, с. 5-19.

[11] Подкопаев С.А., Гаврюшин С.С., Николаева А.С. Анализ процесса нелинейного деформирования гофрированных мембран. Сб. тр. Математическое моделирование и экспериментальная механика деформируемого твердого тела, 2017, вып. 1, с. 31-36.

[12] Подкопаев С.А., Гаврюшин С.С., Николаева А.С., Подкопаева Т.Б. Расчет рабочей характеристики перспективных конструкций микроактюаторов. Сб. тр. Математическое моделирование и экспериментальная механика деформируемого твердого тела, 2017, вып. 1, с. 45-51.

[13] Пономарев С.Д., Бидерман В.Л., Лихарев К.К. и др. Расчеты на прочность в машиностроении. Т.2. Москва, Машгиз, 1958, 975 с.

[14] Феодосьев В. И. Упругие элементы точного приборостроения. Москва, Оборонгиз, 1949, 342 с.

[15] Феодосьев В.И. К расчету хлопающей мембраны. Прикл. математика и механика, 1946, т. 10, № 2, с. 295-306.

[16] Гаврюшин С.С. Численное моделирование процессов нелинейного деформирования тонких упругих оболочек. Математическое моделирование и численные методы, 2014, № 1, с. 115-130.

[17] Belhocine A. Exact analytical solution of boundary value problem in a form of an infinite hypergeometric series. International Journal of Mathematical Sciences and Computing (IJMSC), 2017, vol. 3, № 1, pp. 28-37.

DOI: 10.5815/ijmsc.2017.01.03.

[18] Crisfield M.A. A fast incremental/iterative solution procedure that handles "snapthrought". СотрШ. and Structures, 1981, vol. 13, № 1, pp. 55-62.

[19] Chuma F.M., Mwanga G.G. Stability analysis of equilibrium points of newcastle disease model of village chicken in the presence of wild birds reservoir. International Journal of Mathematical Sciences and Computing (IJMSC), 2019, vol. 5, № 2, pp. 1-18. D0I:10.5815/ijmsc.2019.02.01

[20] Gupta N.K., Venkatesh. Experimental and numerical studies of dynamic axial compression of thin walled spherical shells. Int. J. of Impact engineering, 2004, vol. 30, pp. 1225-1240.

[21] Marguerre, K.: Zur: Theorie der gerkrümmten Platte großer Formänderung. Proceedings of the Fifth International Congress for Applied Mechanics, 1939, pp. 93-101.

[22] Moshizi M.M., Bardsiri A.K. The application of metaheuristic agorithms in automatic software test case generation. International Journal of Mathematical Sciences and Computing (IJMSC), 2015, vol.1, №3, pp. 1-8.

DOI: 10.5815/ijmsc.2015.03.01

[23] Mescall J. Numerical solution of nonlinear equations for shell of revolution. AIAA J., 1966, vol. 4, № 11, pp. 2041-2043.

[24] Reissner E. Alisymmetrical deformations of thin shells of revolution. Proc. of Symp. In Appl. Math., Ameг. Math. Soc., 1950, vol. 3, pp. 27-52.

[25] Riks E. The application of Newton's method to the problem of elastic stability. J. Appl. Mech., 1972, vol. 39, pp. 1060-1065.

Статья поступила в редакцию 14.12.2019

Ссылку на эту статью просим оформлять следующим образом: Подкопаев С.А. Численное моделирование закритического нелинейного деформирования осесимметричных мембран. Математическое моделирование и численные методы, 2020, № 1, с. 64-87.

Подкопаев Сергей Анатольевич — канд. техн. наук, доцент кафедры «Компьютерные системы автоматизации производства» МГТУ им. Н.Э. Баумана. Область научных интересов: численные методы и алгоритмы анализа процессов нелинейного деформирования тонкостенных элементов конструкций машин и приборов, аддитивные технологии и прототипирование, автоматизация технологических процессов и производств. e-mail: [email protected], [email protected]

Numerical simulation the post-buckling nonlinear deformation of axisymmetric membranes

© S.A. Podkopaev Bauman Moscow State Technical University, Moscow, 105005, Russia

The theoretical foundations of nonlinear straining of thin-walled axisymmetric shells are considered. The operational characteristics of the membranes in various switching devices, valves and pressure sensors are presented. The types of non-linear behavior ofpost-buckling behavior of axisymmetric membranes are considered. A mathematical model is presented to describe nonlinear straining of axisymmetric membranes, a discrete continuation by parameter method, and the "changing the subspace of control parameters" technique. Using the hinged spherical shell as an example, a study of post-buckling behavior is performed. A rational mathematical model has been selected to describe nonlinear straining of thin-walled axisymmetric shells. A numerical algorithm for studying the processes of nonlinear straining of multi-parameter systems has been developed and implemented as an author program.

Keywords: nonlinear straining, thin-walled axisymmetric shell, membrane, post-buckling behavior, discrete switching, continuation by parameter, change of the subspace of parameters

REFERENCES

[1] Alfutov N.A. Osnovy rascheta na ustojchivost' uprugih sistem. Biblioteka raschetchika [Fundamentals of calculation on the stability of elastic systems. Calculator library]. Moscow, Mashinostroenie Publ., 1977, 488 p.

[2] Andreeva L.E. Uprugie elementy priborov: uchebnoe posobie [Elastic elements of devices: Tutorial]. Moscow, Mashinostroenie Publ., 1982, 456 p.

[3] Biderman V.L. Mekhanika tonkostennyh konstrukcij. Statika. Biblioteka raschetchika [Mechanics of thin-walled designs. Statics. Calculator library]. Moscow, Mashinostroenie Publ., 1977, 488 p.

[4] Valishvili N.V. Metody rascheta obolochek vrashcheniya na ECVM [Methods for calculating the shells of rotation on electronic digital computer]. Moscow, Mashi-nostroenie Publ., 1976, 278 p.

[5] Volmir A.S. Ustojchivost' deformiruemyh sistem [Resistance of deformable systems]. Moscow, Fizmatgiz Publ., 1967, 984 p.

[6] Gavrushin S.S. Razrabotka metodov rascheta i proektirovaniya uprugih ob-olochechnyh konstrukcij pribornyh ustrojstv: dissertaciya [Development of methods for calculating and designing elastic shell structures of instrument devices: dissertation]. Moscow, 1994, 316 p.

[7] Gavrushin S.S., Baryshnikova O.O., Boriskin O.F. Chislennye metody v dinamike i prochnosti mashin [Numerical methods in dynamics and strength of machines]. Moscow, BMSTU Publ., 2012, 492 p.

[8] Grigolyuk E.I., Lopanitsyn E.A. Konechnye progiby, ustojchivost' i zakriticheskoe povedenie tonkih pologih obolochek [Finite deflections, stability and post-buckling behavior of thin shallow shells]. Moscow, Moscow State University of Mechanical Engineering Publ., 2004, 162 p.

[9] Grigolyuk E.I., Kabanov V.V. Ustojchivost' obolochek [The stability of the shells], Moscow, Nauka Publ., 1978, 359 p.

[10] Il'gamov M.A. Trudy Matematicheskogo centra im. N.I. Lobachevskogo — Proceedings of the Lobachevsky Mathematical center, 2010, vol. 42, pp. 5-19.

[11] Podkopaev S.A., Gavrushin S.S., Nikolaeva A.S. Analiz processa nelinejnogo de-formirovaniya gofrirovannyh membrane [Analysis of the process of nonlinear strain-

ing of corrugated membranes]. Sb. tr. Matematicheskoe modelirovanie i eksperi-mental'naya mehanika deformiruemogo tverdog [Collection of works. Mathematical modeling and experimental mechanics of a deformable solid], 2017, iss. 1, pp. 31-36.

[12] Podkopaev S.A., Gavrushin S.S., Nikolaeva A.S., Podkopaeva T.B. Raschet rabochej harakteristiki perspektivnyh konstrukcij mikroaktyuatorov [Calculation of the working characteristics of promising designs microactuators]. Sb. tr. Ma-tematicheskoe modelirovanie i eksperimental'naya mehanika deformiruemogo tver-dog [Collection of works. Mathematical modeling and experimental mechanics of a deformable solid], 2017, iss. 1, pp. 45-51.

[13] Ponomarev S.D., Biderman V.L., Likharev K.K. et al. Raschety na prochnost' v mashinostroenii. T.2. [Strength calculations in mechanical Engineering]. Moscow, Mashgiz Publ., 1958, 975 p.

[14] Feodosyev V.I. Uprugie elementy tochnogo priborostroeniya [Elastic elements of precision instrument engineering]. Moscow, Oborongiz Publ., 1949, 342 p.

[15] Feodosyev V.I. Prikl. matematika i mekhanika — Applied Mathematics and Mechanics, 1946, vol. 10, no. 2, pp. 295-306.

[16] Gavriushin S.S. Matematicheskoe modelirovanie i chislennye metody — Mathematical modeling and Computational Methods, 2014, no. 1, pp. 115-130.

[17] Belhocine A. Exact analytical solution of boundary value problem in a form of an infinite hypergeometric series. International Journal of Mathematical Sciences and Computing (IJMSC), 2017, vol. 3, no. 1, pp. 28-37. DOI: 10.5815/ijmsc.2017.01.03

[18] Crisfield M.A. A fast incremental/iterative solution procedure that handles "snapthrought". Comput. and Structures, 1981, vol. 13, no. 1, pp. 55-62.

[19] Chuma F.M., Mwanga G.G. Stability analysis of equilibrium points of newcastle disease model of village chicken in the presence of wild birds reservoir. International Journal of Mathematical Sciences and Computing (IJMSC), 2019, vol. 5, no. 2, pp. 1-18. DOI:10.5815/ijmsc.2019.02.01

[20] Gupta N.K., Venkatesh. Experimental and numerical studies of dynamic axial compression of thin walled spherical shells. Int. J. of Impact engineering, 2004, vol. 30, pp. 1225-1240.

[21] Marguerre, K.: Zur: Theorie der gerkrümmten Platte großer Formänderung.

Proceedings of the Fifth International Congress for Applied Mechanics, 1939, pp. 93-101.

[22] Moshizi M.M., Bardsiri A.K. The application of metaheuristic agorithms in automatic software test case generation, International Journal of Mathematical Sciences and Computing (IJMSC), 2015, vol.1, no. 3, pp. 1-8. DOI: 10.5815/ijmsc.2015.03.01

[23] Mescall J. Numerical solution of nonlinear equations for shell of revolution. AIAA J., 1966, vol. 4, no. 11, pp. 2041-2043.

[24] Reissner E. Alisymmetrical deformations of thin shells of revolution. Proc. of Symp. In Appl. Math., Ameг. Math. Soc., 1950, vol. 3, pp. 27-52.

[25] Riks E. The application of Newton's method to the problem of elastic stability. J. Appl. Mech, 1972, vol. 39, pp. 1060-1065.

Crisfield M.A. A fast incremental/iterative solution procedure that handles "snapthrought". Comput. and Structures, 1981, vol. 13, № 1, pp. 55-62.

Podkopaev S.A., Cand. Sc. (Eng.), Assoc. Professor, Department of Computer systems for production automation, Bauman Moscow State Technical University. Research interests: numerical methods and algorithms for analyzing the processes of nonlinear deformation of thin-walled structural elements of machines and equipment, additive technologies and prototyping, automation of technological processes and production. e-mail: [email protected], [email protected]