Расчет нелинейного напряженно-деформированного состояния непологих оболочек вращения под действием ветровой нагрузки Текст научной статьи по специальности «Физика»

УДК 539.37

РАСЧЕТ НЕЛИНЕЙНОГО НАПРЯЖЕННО-ДЕФОРМИРОВАННОГО СОСТОЯНИЯ НЕПОЛОГИХ ОБОЛОЧЕК ВРАЩЕНИЯ ПОД ДЕЙСТВИЕМ ВЕТРОВОЙ НАГРУЗКИ

М.С. ГАНЕЕВА, В.Е. МОИСЕЕВА, З.В.СКВОРЦОВА

Институт механики и машиностроения КазНЦ РАН

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

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

На предприятиях топливно-энергетического комплекса широко используются резервуары разного рода для хранения жидких и газообразных продуктов, а также для технологических целей. В большинстве аппаратов объем ограничен цилиндрической обечайкой с крышками сферической, эллиптической или конической формы. Емкостями такого типа являются корпуса газгольдеров, резервуары для топлива и нефтепродуктов, баки для осветления воды в схемах химической водоочистки ТЭС, резервуары для хранения сжиженных газов и др. В виде цилиндрических оболочек со сферическим куполом выполняются защитные оболочки реакторов АЭС. В процессе эксплуатации купольные конструкции подвергаются различным нагрузкам: давлению внутренней среды, перепаду температур, ветровой, снеговой нагрузке.

Данная работа посвящена исследованию воздействия нагрузки типа ветровой на купол крупного резервуара, выполненный в виде оболочки вращения с центральным отверстием. В реальном изделии отверстие может быть закрыто люком или соединено с вентиляционным патрубком. Расчет представляет трудности в случае интенсивного нагружения с достижением предельных нагрузок потери устойчивости. В работе дается обобщение алгоритма авторов [1] для расчета геометрически и физически нелинейного напряженно-деформированного состояния непологих оболочек вращения при неосесимметричном нагружении. Основу методики составляет алгоритм продолжения численного решения по некоторому интегральному параметру, который определяется эпюрой полуволны на амплитуде одной из основных гармоник прогиба в тригонометрическом ряду по окружной координате. Предложенный алгоритм позволяет получать при численных расчетах предельные точки на характерных кривых нагружения (критические значения параметра неосесимметричной нагрузки). Представлены результаты расчетов для полусферической оболочки с центральным отверстием с жестко заделанными краями под действием неосесимметричного нагружения типа ветровой нагрузки. Задаваемая нагрузка моделирует давление ветра на купол вертикального резервуара.

1. Рассматривается оболочка вращения (рис. 1), замкнутая в окружном направлении [2, 3]. За координатные линии приняты меридианы 5 , параллели ^ и

7 ТТ «1 < 5 < ,0 <ф< 2п,

внешняя нормаль 5 к поверхности приведения. При этом 1 Н ' ^ '

-Н1 < 7 < к2, й(«) = Н1 + Н2 - А1(«), к2(«) -

1 * 2» V/ 1 2 толщина оболочки; ^ " 2 главные кривизны;

0 £ [0 0 ] —

г - радиус параллели; 1 1 угол между осью вращения оболочки и

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

xi (ф)? i 1,3 и температурное поле T(s'фz) при начальной температуре

T° — const. Материал оболочки изотропный с характеристиками, зависящими от температуры. Задача решается на основе уравнений теории Кирхгофа-Лява в геометрически нелинейной постановке при умеренных поворотах. Влияние температурного поля учитывается по гипотезе Дюгамеля-Неймана. Допускается работа материала за пределом упругости, связь между напряжениями и деформациями описывается уравнениями малых упругопластических деформаций без учета разгрузки. Основные соотношения получены в недеформированных координатных линиях. Далее используются соотношения и обозначения работы [2].

Рис. 1. Оболочка вращения

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

Xf (s, ф) — (s)cos кф, i —1»3;

к—0

X 21 (s, ф) — £ X£k (s)sin кф,

к—1

n

Т - То = ^ Тк (I )соз кф; 0 <ф<п. к=0

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

У = (Т1н1,Т1н2, ен + -М", Мп, П, ъ, »!)' г дф

в тригонометрические ряды по окружной координате [4]:

Ь Ь

™ =Х^(*)со«кф, V2 = £V2,k(«)«1Икф, К

к=0 к=1 (1)

В силу нелинейности задачи Ь > п, и значение Ь устанавливается численным экспериментом.

Вводится вектор разрешающих функций для амплитуд рядов (1):

2к = (Т11,к ,Т12,к, е1,к, М11,к, %к, v2,k, ™к, ®1,к )', к = 0 Ь. (2)

Задача сводится к интегрированию ряда систем для нелинейных обыкновенных дифференциальных уравнений:

Л2к

= Нк («) г к + Як (,Тк) + Гк (s,у) + Фк (у), к = 0, Ь ж (3)

при соответствующих граничных условиях:

Акгк = ак при s = вкгк = ьк при s = ^, к = 0Ь. (4)

Здесь Нк (s)- матрица коэффициентов (8 х 8); Як - вектор (8 х 1) амплитуд

нагрузочных и температурных членов; Ак,Вк - матрицы (4 х 8); ак,Ьк - векторы (4 х 1). При этом коэффициенты Фурье разложений связывающих системы (3)

нелинейных членов к («, ), к («, ) в ряды вида (1) подсчитываются численно по формулам Бесселя [5].

Рассмотрим случай, когда в качестве ведущего параметра задачи выбран параметр нагружения (параметр нагрузки или температуры). При данном значении ведущего параметра, номер которого в дальнейшем опущен, нелинейная краевая задача (3), (4) решается методом последовательных приближений в сочетании с методом ортогональной прогонки [6]. В качестве метода последовательных приближений реализованы метод общей итерации

<<7 (т+1)

=нк (+1)+Ек («, х»к, тк)+л(т >,

¿« к к (5)

л(кт) = Гк («,¥(т ))+Фк («,¥(т)), к=01 , т > 0

и метод линеаризации с подсчетом матрицы Якоби только от геометрически

Г 7

нелинейных членов к по вектору к :

<<7 (т+1)

к = щ(т) + щ(т) + щ (т) <« "Щ1'к +Щ2'к +Щ3'к , (6)

щ(т ) щ1,к

щ («)+«Гк<«, 70" ),А.7? >)

7(т+1) 7к ,

<7к

= Ек(«,Х*к,Тк) + Гк(«,70т),л ,7[т)),

(т) (т)

щ^ = - <Гк(«,7,<<7,л ,7, ) 7т) + фк(«,70т),К 7(т)),

к = 0, Ь, т > 0.

За ¥(0) принимается решение, полученное с заданной точностью на

у(0) = 0 П 7кт+1), к = 0Ь

предыдущем шаге, или 1 — ". По полученным к подсчитываются

промежуточные векторы ~(т+1), затем последующее приближение

¥(т+1) = ¥(т) + т(~(т+1) - ¥(т)), 0 < т < 1 (7)

с использованием коэффициента релаксации т и амплитуды нелинейных членов

Г (« ¥(т+1)) Ф (« ¥(т+1)) ^^ ', КК ' ' по формулам Бесселя [5]. Процесс продолжается до

достижения заданной точности 5 по какому-либо критерию.

Процесс последовательных приближений (5)-(7), когда ведущим параметром

выбран параметр нагружения, обеспечивает получение решения только на монотонно

возрастающем участке характерной кривой нагружение-прогиб. Предлагается

использовать идею о введении в качестве ведущего параметра некоторой

интегральной величины, ранее разработанную в работах [7-9] для осесимметричных

задач, при решении уравнений неосесимметричной задачи (3), (4).

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

к , полагая ™к амплитудой одной из основных гармоник в ряду (1). Введем интегральную амплитуду [1]:

Ук = } ^, Ук (Sl) = 0, Ук (SN) = Ск.

41 (8)

Неизвестные 0 , Ук удовлетворяют уравнениям: — = 0; —к = >~к.

ds ds (9)

Введем вектор неизвестных (10 х 1):

Вк = (ТС11,к ,Т12,к, О1,к, ММ11,к, С1,к, С2,к, Мк, С1,к, О,Ук)', соответствующими уравнениями для которого будут (3) при к = к и (9):

^ = ННк (s) Як + (s,Tk)+Г к (s,у) + Ф к (s,у), ^ (10)

где Нк (s)- матрица коэффициентов (10 х 10), остальные величины - векторы (10 х 1). Граничные условия (4), (8) запишутся в виде:

АкЯк = ак при s = ВкЯк = Ск при S = ^, (11)

где Ак,Вк - матрицы (5 х 10), ак,Ьк - векторы (5 х 1).

В краевой задаче (10), (11) ведущим параметром процесса нагружения могут

выступать величины 0 или Ск. В случае нарастания в процессе деформирования

всей эпюры амплитуды Wk монотонное увеличение параметра Ск обеспечит

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

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

выделить ту часть Si < < ^ , на которой будет наблюдаться монотонное возрастание

величины интегральной амплитуды прогиба °к . В таком случае алгоритм на основе (8), (9) допускает обобщение. Пусть меридиан оболочки разбит на три части следующим образом:

Ук = а11М к^ йУ^ =

Sl < s < si, Sl ^ а1Мк,

(12

а)

ук (s1) = 0, Ук(ч) = с к (1);

Si < s < s¡,

Ук = а2 / + ск(1), йУСк =

ds

= а2 Мк >

б)

в)

Ук() = ск(1)' Ук) = ск(2);

5

Ук = «3 / + С к (2), ¡У~к ~ '] * 5 ^ ^ , Ч -¡к = '

(12

Ук (Я] ) = Ск (2)' ук (Ы ) = Ск (3).

Некоторые случаи применения алгоритма (12а)-(12в) в зависимости от значений коэффициентов а1' «2' «3 показаны в табл. 1.

Таблица 1

Варианты применения алгоритма

«1 «2 «3 ведущий параметр

1 1 1 Ск(3) = Ск

1 0 0 Ск (3) = Ск(2) = Ск(1)

0 1 0 ск(3) = Ск(2)

0 0 1 Ск(3)

Кратко опишем алгоритм численного решения нелинейной краевой задачи (3)-

(4). ~

1. Определим в качестве ведущего параметра задачи значение Ск (3)

интегральной амплитуды Ук на правом краю интервала 5 = , с учетом данных табл. 1. В алгоритме расчетов для системы (10) используются процессы (5) или (6).

Далее для краткости обозначим к(3) .

2. Решаем краевую задачу (10), (11) без учета нелинейных членов, подбираем

начальное значение параметра С и его шаг изменения . Используя С, решаем линейную краевую задачу (10), (11), получаем вектор ' , в который входит первое

0(1) с

приближение параметра нагрузки 1 на шаге 1.

0(1) У(1) Гк (5 У(1))

3. Передаем ^ в краевую задачу (4)-(7), получаем 1 , кУ ' 1 , К1)

Фк («УГ), к = о, ь.

Л ТТ " Гк (У(1) ), Фк (У(1)) „„

4. Передаем нелинейные члены Л 4,1 «V'1 / в задачу (10), (11), получаем вектор , в том числе 1 .

5. Передаем 01 в краевую задачу (4)-(7), получаем У1 , Г к (^^ ),

Фк(5,У1(2)), к = 0,ь.

6. Повторяем пп. 4, 5 с нарастанием номера последовательных приближений

р У (р) = У

до достижения заданной точности за итераций. Получаем решение 1 1 для

C Q( Р) — Q.

параметров интегральной амплитуды 1 и нагрузки 1 1.

C — C + AC

7. Переходим к следующему шагу параметра 2 — 1 , для которого за

Y (0) — Y

нулевое приближение принимается 2 1.

8. Повторяем пп. 4, 5, 7 до получения Yi' Qi' i — 11 и других характеристик напряженно-деформированного состояния задачи.

Алгоритм реализован в виде программы на языке ФОРТРАН для ПЭВМ.

2. Рассмотрим деформирование сферического купола с центральным отверстием под действием нагрузки типа ветровой. Обозначим R и h - радиус сферы

Rlh —100 0i — 0,2, 0N —я/2 тл й и толщину; ' , 1 ' ' N 1 . Края оболочки жестко заделаны:

V1 — V2 — w — »1 — 0 при s — s1, s — sn .

Материал оболочки следует закону линейного упрочнения [10] с коэффициентом упрочнения ^ — 0,9 ; v — 0,3 - коэффициент Пуассона,

_3

slE — 7,459'10 - предел текучести, E - модуль Юнга, — ( 2 2 _ 3 2 )0,5

CTi — (СТ11 + °22 _ СТ11СТ22 + °12) - интенсивность напряжений. На оболочку действует внешнее нормальное давление

X3 — Q(1 + cos ф + 0'5cos2ф) (13)

В формуле (13) параметр нагрузки Q представляет собой среднее давление, которому пропорционален весь груз на поверхность оболочки. При оссесимметричном нагружении

X3н = Q. (14)

Численные расчеты проведены на основе алгоритма, описанного в п. 1, при

к — 0, L — 8, N — 400, в —10_4 , ч

с выбором параметра процесса (12а)-(12в) при

a\ —1' a2 — аЪ — 0' si — s65. Результаты вычислений приведены на рис. 2-5. Здесь приняты сокращения: Л - линейное решение; Г - геометрически нелинейное решение; Ф - физически нелинейное решение; ГФ- геометрически и физически нелинейное решение. Сплошные линии отражают неосесимметричное нагружение (13), штриховые - осесимметричное нагружение (14).

По зависимостям параметра

Рис. 2. Зависимости параметра нагружения Q

, нагружения ^ от максимального

от максимального прогиба rj w

прогиба w, полученным при решении задачи в Л-, Г-, Ф-, ГФ- постановках (рис. 2), видно, что уже при малых

- \wlh»0,5 -прогибах наблюдается

количественное различие в

результатах. Здесь же проявляется и

качественное влияние совместного

учета ГФ-нелинейностей: ГФ-решение

выявляет предельные нагрузки как при

осесимметричном, так и

неосесимметричном нагружении. На

ГФ-кривой при неосесимметричном

нагружении отмечены номерами пять

характерных точек, которым

соответствуют состояния: 1) линейная

Q/E = - 2,79 ■ 10_5; 3)

_5 max a i = a S

область,Q' = _ 2,45 ■10 ; 2) достигнут предел текучести материала z -5 ч

" достигнуто верхнее критическое значение параметра

неосесимметричной нагрузки QIE = _ 3,37 ■10 ; 4) ниспадающая ветка, QIE =

_ 2 39 10_5; QlE _ 210 10_5

' , 5) ниспадающая ветка, = • Верхнее критическое значение

параметра осесимметричной нагрузки составляет QIE = _10,72 ■10 (точка 6), что примерно в 3 раза выше соответствующего значения в неосесимметричном случае. На рис. 3 показаны зависимости интенсивности напряжений от параметра нагрузки

на наиболее нагруженном меридиане ^ = 0. Отмечено безразмерное значение предела

a s = 103 а S/e .

текучести материала ^ — М • Видно, что в ГФ-решении кривые зависимостей

стг(0) изменяются немонотонно: после достижения предела текучести рост

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

Рис. 3. Зависимости интенсивности напряжений от параметра нагрузки

На рис. 4, 5 представлены эпюры по меридиану ^ = 0, полученные в ГФ-решении для пяти состояний, отмеченных номерами на рис. 2. На рис. 4, а показан

прогиб Отмечен интервал на котором определяется ведущий параметр

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

Wо' тригонометрического ряда (1). Эти графики показывают достаточность

заданного количества членов Ь в ряду. На рис. 5 даны эпюры интенсивности

напряжений ст 1 на внутренней г = —0и внешней г = лицевых поверхностях. Эпюры носят сложный характер. С нарастанием процесса деформирования превышается предел текучести материала.

Рис. 4. Эпюры прогиба ™ и амплитуды ™0 ' ™1 ' ™ 8 тригонометрического ряда (1)

Рис. 5. Эпюры интенсивности напряжений

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

Summary

Effects of wind-like loading on the dome of large reservoir designed as a shell of revolution with central hole are investigated, in case of intensive loading when ultimate loads of instability are reached. The algorithm of computation of geometrically and material nonlinear stress-strained state allowing to receive critical parameter values of non-axisymmetric loading is proposed. Results of computations for hemispherical dome are

submitted and analyzed.

Key words: stress-strained state, shell of revolution, instability, ultimate loads.

Литература

1. Ганеева М.С., Моисеева В.Е. Методика расчета больших прогибов непологих упругопластических оболочек вращения при неосесимметричном нагружении // Известия вузов. Авиационная техника. 2007. № 4. С. 3-7.

2. Ганеева М.С. Термосиловая задача в геометрически и физически нелинейной теории нетонких и тонких оболочек / КФТИ КФАН СССР. Казань, 1985. 126 с. Деп. в ВИНИТИ 24.06.85. № 4459-85Деп.

3. Ганеева М.С. Прочность и устойчивость оболочек вращения. М.: Наука, 1992. 161 с.

4. Ганеева М.С., Косолапова Л.А., Моисеева В.Е. Расчет на прочность гибких упругопластических оболочек вращения при неосесимметричном термосиловом нагружении // Актуальные проблемы механики оболочек. Труды Междунар. конф. Казань: УНИПРЕСС. 1998. С. 35-42.

5. Ланцош К. Практические методы прикладного анализа: Справочное руководство. М.: Физматгиз, 1961. 524 с.

6. Годунов С.К. О численном решении краевых задач для систем линейных обыкновенных дифференциальных уравнений // Успехи математич. Наук. 1961. 16. № 3. С. 171-174.

7. Ганеева М.С., Алексеева О.В. Об одном алгоритме численного решения геометрически нелинейных осесимметричных задач непологих оболочек вращения // Исследования по теории оболочек: Труды семинара. Казан. физ.-техн. ин-т КФ АН СССР, 1976. № 7. С. 120-127.

8. Ганеева М.С. Нелинейный осесимметричный изгиб непологой оболочки вращения средней толщины // Прочность и устойчивость оболочек: Труды семинара. Казан. физ.-техн. ин-т КФ АН СССР, 1980, № 13. С. 29-41.

9. Ганеева М.С., Моисеева В.Е. Нелинейный изгиб нетонких составных оболочек вращения из термочувствительного упругопластического материала // Исследования по теории оболочек: Труды семинара. Казан. физ.-техн. ин-т КазНЦ АН СССР, 1990. № 25. С. 4-20.

10. Ильюшин А.А. Пластичность. Ч.1. Упругопластические деформации. М.-Л.: Гостехтеориздат, 1948. 376 с.

Поступила в редакцию 10 сентября 2008 г.

Ганеева Музайна Саитгареевна - д-р физ.-мат. наук, старший научный сотрудник, главный научный сотрудник Учреждения Российской академии наук Института механики и машиностроения Казанского научного центра РАН (ИММ КазНЦ РАН).

Моисеева Валерия Евгеньевна - канд. физ.-мат. наук, старший научный сотрудник Учреждения Российской академии наук Института механики и машиностроения Казанского научного центра РАН (ИММ КазНЦ РАН).

Скворцова Зара Владимировна - канд. физ.-мат. наук, ученый секретарь Учреждения Российской академии наук Института механики и машиностроения Казанского научного центра РАН (ИММ КазНЦ РАН). Тел. 292-51-62. E-mail zara-skvortsova@yandex.ru.