Численное моделирование динамики осесимметричного раскрытия мягкой оболочки поле силы тяжести Текст научной статьи по специальности «Механика и машиностроение»

№ 10 2006 РАСЧЕТ И КОНСТРУИРОВАНИЕ МАШИН

539.3

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ДИНАМИКИ ОСЕСИММЕТРИЧНОГО РАСКРЫТИЯ МЯГКОЙ ОБОЛОЧКИ

В НОЛЕ СИЛЫ ТЯЖЕСТИ

Канд. техн. наук, дои. П.Г. РУСАНОВ

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

There is wrote the simplest manner of numerical analysis for dynamic process expansion soft flexible shell in gravity space during gas filling under pressure by mean the discretization technique solid bodies methods. For axes-symmetrical case of test there is obtained the numerical information about shell static condition, when pressure is absent, and also about the next changes of meridian form, thickness and tensions on account of pressure rising.

Теоретический анализ скоротечного динамического процесса раскрытия замкнутой мягкой оболочки при резком повышении внутреннего давления представляет собой нелинейную задачу механики для системы с бесконечным числом степеней свободы. Предлагается решение осесымметричгюго варианта задачи на основе дискретной расчетной схемы метода твердых тел (МТТ) [1—12], который, по сравнению с MIO, в задачах такого класса имеет более наглядные расчетные схемы, простой вывод уравнений состояния, технологичные алгоритмы численного решения и лучшую точность результатов.

Постановка задачи. Симметричная замкнутая область контейнера (рис. 1) ограничена однородной тонкой оболочкой из эластика (коэффициент Пуассона v =0,5) толщиной h массой М0 + 2ЛМ с жестким подвижным дном массой М. Оболочка подвешена на неподвижной вертикальной трубе — штуцере и наполняется газом. В ненагруженном состоянии ее срединная поверхность имеет форму круглого цилиндра радиусом R и длиной образующей L0+ AL. Внешние радиусы дна и трубы одинаковы и равны R - 0,51г(У До начала раскрытия (время t < 0) конструкция неподвижна в поле силы тяжести при нулевом избыточном давлении газа р внутри оболочки. Главные допущения: безмоментное напряженное состояние материала оболочки с постоянной удельной плотностью, пренебрежение массой газа внутри оболочки, тепловыми процессами и движением воздуха снаружи контейнера. Цель — в рамках наглядных физических представлений изучить переходные процессы изменения формы меридиана, толщины и напряженного состояния оболочки при известном законе роста давления р = р(г).

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

№ 10

2006

X

Ш О 1

Р

2

£

Рис. 1. / — оболочка, 2 — дно, 3 — труба

ностей (плоскость и конус), образованных движением нормали к срединной поверхности оболочки по контурам «вживленной» в нее сетки. В кеиагруженпом состоянии оболочки линии сетки расположены в равноотстоящих п + 2 горизонтальных и п0 аксиальных плоскостях. Поэтому в этом состоянии все (п + 2) х 2п0 элементов одинаковы, имея массу т = 0,5 Мс/(пп0) и характерные размеры 1() х 2/?(3 х 1г0.1(} = Ь0/(п + 1), ¡3 = 0,5п/по.

Г # с н I

НИ 1

Рис. 2

Упростим граничные условия искомой математической модели, согласовав ранее назначенные величины АЬ, М4о с шагом сетки разбиения, а именно, положим М - 1() и ДМ = 2тп(У Тогда ДМ — масса, а ДЬ — ширина одного пояса ненагруженной оболочки. По той же причине считаем, что крайние пояса с номерами к - 0 и к = п + 1 на половину своей ширины 0,5/„ заходят на боковые цилиндрические поверхности дна и трубы и жес-

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

№ 10 2006

тко связаны с ними, а центры масс крайних поясов всегда лежат в соответствующих торцовых плоскостях тел закрепления (г0 = 0, z — аппликаты торцов трубы и дна).

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

О движении каждого элемента будем судить по изменениям положения его собственного репера, т.е. главных центральных осей инерции массы элемента. С учетом осесимметричности формы оболочки обобщенными координатами системы «дно, реперы элементов оболочки» могут служить z — аппликата торца дна и тройки однотипных геометрических параметров, задающих положения реперов для элементов внутренних поясов с номерами к ~ 1,2,..., п: 1) хк — расстояние от центра масс элемента (точки Ск) до оси Oz; 2) zk, — аппликата центра масс в неподвижной системе координат; 3) ак — угол наклона оси инерции СкС, элемента к плоскости параллели.

Чтобы напрямую следить за формой оболочки, вместо координаты хк будем использовать rk = ОкЛк — радиус параллели срединной поверхности оболочки в плоскости, проходящей через Ск — центр масс элемента (Ок — центр кривизны параллели, Ак— точка параллели, zAk = zk, рис. 2, б). Упростим связь между хк и гА„ пренебрегая на данном этапе размерами поперечного сечения элемента в аксиальной плоскости, проходящей через ось Oz, Тогда, приближенно, хк ~ sinp гк/р, причем точность выполнения данного соотношения повышается с ростом п — числа внутренних поясов. В результате получаем возможность указывать положение дискретных точек меридиана Ак по значениям координат точек Ск, рассчитываемых на основании уравнений динамики.

Уравнения движения. В рассматриваемых условиях у внутренних поясов репер каждого элемента имеет три степени свободы. Однако в целях упрощения применим лишь уравнения движения для координат rk, zk, полагая, в силу допущения о безмомент-ном состоянии оболочки, движения реперов по координатам ак быстрыми, а по координатам гк, zk медленными. Вместо ак — угла наклона оси репера СкС, будем следи ть за Qk — углом направления касательной к меридиану оболочки в его дискретных точках Ак, оценивая tgQk по отношению конечных разностей кординат rk, zk

tg6A, = Zk+l ~ Vl , (* = 1,2,..., п). (1)

Гк+1 Гк-1

Заметим, что ак —> при бесконечном уменьшении параметров р и 1(). В точках А0, Ап+{ значения углов Qk доопределяем формулами 260=0,57i + 0р 20/г+1 = 0,571 + 0/Г

Поскольку все элементы одного пояса имеют одинаковые значения координат rk, zk, упростим задачу, сократив число элементов в каждом поясе до двух. Тогда 2л0 = 2, 2(3 - тх, масса элемента т - 0,5М/п, хк =2/пгк.

Чередуя номера поясов к = 1,2и, составим по два скалярных уравнения динамики на основе теоремы о движении центра масс, примененной к массе лишь одного из элементов пояса. В них применим упрощенные аналитические выражения для проекций главных векторов внешних распределенных сил, получаемые с использованием аппроксимации линии меридиана срединной поверхности—ломаной А0, А1 , ... Ак, Ап>Ап+[ . Для элемента к-ro пояса выражения проекций будут следующими: на радиальную ось, проходящую через центр масс элемента,

Nk.r =2(аА',ш2 Гк,2 \2COS0U-aAv»l Гк,\ hk,\ COS0u)> Tk,r = % I lk,\ + hk,2 1 k,2 ) ~ 0T

меридиональных и окружных сил в оболочке; Pk r = p(t) Skr — от сил избыточного давления газа; Qk r = -2\xSkrr к — от сил вязкого сопротивления внешнего газа;

№10

2006

на вертикальную ось г

^кс ~ П(®кт2 Гк7 ^к 2 2 ~ °кт\ Гк 1 ¡1к 1 ^ — ОТ МерИДИОНаЛЬНЫХ СИЛ;

Рк 7 = рЬ) Зк7 — от сил избыточного давления газа; Ок ,= | | г к — от сил вязкого сопротивления внешнего газа; Ск 7= — от сил тяжести.

Здесь коэффициент 2 в выражениях для <2А, г получен в результате аналитического интегрирования проекций элементарных распределенных сил по дуге параллели срединной поверхности в пределах границ элемента; дополнительные нижние индексы 1 и 2 указывают на соседство границы или участка элемента, соответственно, с нижним или верхним поясами; =0,5тг (г2к 2- г2к,), = (0,5гк } + гк + 0,5гк + Нк — площади проекций срединной поверхности элемента деформированной оболочки на две плоскости, перпендикулярные, соответственно, оси Ог. и радиальной оси, проходящей через центр масс элемента;

2Нк = " ' 2ги = гк + г,_,; 2гкЛ = г,+! + г,;

26 к.\ = + 1 ' 2 = + 0*+1 ^

Кг Уо = 0,5/г0/07?; (2)

21и = V~ гА-1 )2 + (г* - VI )2; 2/„ = V [г,+1 - г, )2 + - г, )2, где 21 к ] , 2/а.2 — длины отрезков ломаной линии, аппроксимирующей меридиан, Нк,, 1гк2 — текущие толщины к-го пояса оболочки по его верхней и нижней границам; окт, ак1, ек т, еА. 1 — осредненные по толщине оболочки напряжения и линейные деформации материала на границах элемента в меридиональных и окружных направлениях, екп — линейная деформация оболочки в нормальном к ее поверхности направлении

= + Ч^«.! ^к,12 = еАЧГ2 + ЧЧш2 + (3)

где Я, = £'/(1-V2); |я > 0 — постоянный коэффициент вязкого сопротивления.

Относительные деформации материала оболочки гш, е,, е на границах элемента &-го пояса оцениваем по изменениям характерных размеров элемента

£к,п=ги/К-Ъ Чм=К\/1га-Ъ

(4)

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

-тг = Рг+Ыг+Т,+{2Г; = Л + + (5)

я

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

(М + 2ш)гд = (М + 2т)g + + - ¿д, (6)

и принятый закон изменения давления р от 0 до ртах на отрезке времени [0, 7], например,

№ 10 2006

= Рт.м (1 - ь = согш > 0. (7)

где Т~ намеченное время исследования процесса заполнения оболочки. Применив обозначения:

для безразмерных переменных г* г* Л* /*, /?*, а*, т (без указания нижних индексов для п z, /. /г, о ) — Г = /?/•*, г - 7„г\ /г = Л„Л*, / = Ь/, р = , о = Я,о*, I = 7т (0 < т < I);

для безразмерных констант — Л = — , С'' = п , С' = ^

Л ' т ' 2т Е

С" - rcW' с» - 7Z/-;.AV' 7 г „

2/л ' '' шЯ ' ; /„7 ' /,_ ' С - " ш ' 2/» '

<7?тУ лх/н E.RhT' ли/?2'/'

К'' =, л. =---, /С- =-; приведем систему ( 1 )—(7) к безразмер-

(M + 2in)L ' (M + 2m) Lu M+2m

ному виду (символ * y безразмерных параметров опущен, дифференцирование переменных — по безразмерному времени т)

= с;"(>ь/:,. + %А*)Мт) - с>; |>u/u + плЦгк +

+ С"'Г2.АЛ?-тгЛ + V(e,2.i + )lCOs0'Jt -+ v(e„.* +c„,,)]eosOu -

-C'A A JtKu + V(e«,U + V<X,2.* +!;„2,jJ-

z* = - - Û)P(T) - a - r;k

+ + + B„2,,)].vm02Jt -

" с"ЧЛД£...и +v(s,i,t + £«i.t)]'™0u;

гл = ^ + + K?rJttJ,[emlJI + v(E|Ul + c„LJ] - K 'r'ntA ;

tg0t = —/?(т)= 1 схр(-///'т);

2/,., =

lk

2/, 2 - - Г, )2 / л2 + (г,,, - z,)2.

Остальные соотношения из (2)—(4) пригодны и для безразмерной формы, если в них условно принять /? = 1г() ~ 1,1() - \/(п + 1).

Начальные условия. По условиям задачи до начала заполнения оболочки газом скоростные компоненты вектора состояния разрешающей системы уравнений (1)—(7) Х(г) = Равны нулю:

гк =0; гк = 0;г =0; {к - 1, /?). Однако пока остается неизвестной статичная форма срединной поверхности оболочки в однородном поле сил тяжести. Прямой путь расчета значений координат , гк, гк (к = 1, /?) дискретных точек Ак меридиана покоящейся оболочки — решить существенно нелинейные, алгебраические уравнения статики, т.е. частный случай системы (1)—(7) при р = 0 и нулевых скоростях и ускорениях.

№10 . 2006

Чтобы не усложнять алгоритм компьютерного решения, для расчета исходной статичной формы оболочки в поле силы тяжести нами применен динамический подход, учитывающий свойства механических систем с диссипацией энергии. Для этого алгоритм численного решения задачи предусматривает предварительный этап, где проводится расчет статичной формы меридиана срединной поверхности оболочки в поле силы тяжести на основании уже имеющейся системы уравнений (1)—(7), но при р = 0 и переменном ускорении изменяющемся в условном времени скачком или монотонно от нуля до номинального значения #=10 м/с2, с использованием тривиальных начальных условий для вектора состояния, соответствующих покою ненагруженной оболочки: гк - R , Ч = * h > zr =(/i + 1 )/„, hk - h() , гтк- Elk= гпк=0, (к = 1, n). Расчет на первом этапе прекращался при снижении норм скоростей и ускорений вследствие учета сил вязкого сопротивления до заданных уровней точности при установившемся номинальном значении g. Время, затраченное на расчет статичной формы меридиана, на последующих этапах не учитывалось. Полученные на первом этапе численным интегрированием системы уравнений значения координат дискретных точек меридиана восполняют для последующего этапа недостающие компоненты вектора состояния в начальный момент заполнения оболочки газом г = 0 с. С качественных позиций данный способ расчета статического состояния конструкции можно отнести к числу методов последовательного нагружения. На втором этапе решаем ту же задачу Коши, но с переменным давлением газа.

Результаты расчетов. Численный анализ решений выполнен для следующих исходных данных: п = 100, М = 10 т = 10 кг, L0 = R = 200ho= 1 м, g = 10 м/с2; ¡1 = J00 Н-с/м\ Е{ = 1 МПа, p(t) = 1 - ехр(-0,5?) КПа, 7= 10 с.

Расчет на первом этапе формы неподвижно висящей «тяжелой» оболочки вращения при р = 0, показал, что, помимо удлинения образующей и перемещения вниз дна на 6,94 мм, оболочка втянулась по всей высоте внутрь, образуя «шейку» в верхней части (рис. 3. поз. 0). Наибольшее радиальное перемещение срединной поверхности внутрь w = -6,32 мм имеет место на уровне 25-го пояса.

Формы меридиана оболочки w = w(zt !)

-0,1 0,0 0,1 0,2 0,3 0,4

0,0

-0,2

-0,4

-0,6

-0,8

-1,0

-1,2

№ 10 2006

На втором этапе рассчитаны формы меридиана оболочки при ее заполнении газом для дискретных моментов времени t = 0,05; 0,1; 2; 4; 10 с (рис. 3, поз. 1—5). Для приня того уровня сил вязкого сопротивления с ростом давления до стационарного значения р - 1 КПа переходный процесс преобразования формы конструкции носи т квазиапс-риодический характер. При этом наблюдается наибольший прогиб линии меридиана в средней части в радиальном направлении (выпучивание оболочки), сопровождаемый сначала подъемом, а затем опусканием дна. В конечный момент расчета / = Т наибольший прогиб оболочки в радиальном направлении наружу составил w = 36,5 см, а дно опустилось дополнительно вниз на 14,2 см. Отметим, что в расчетах с более низким значением коэффициента сопротивления (ц = 10 Н-с/м3) движения всей линии меридиана оболочки и дна визуально воспринимаются по графической информации, выводимой на экран в процессе счета, как затухающие колебания. Практическое отсутствие волн па расчетных изображениях линии меридиана подтверждает правомочность применения метода разделения быстрых и медленных движений при формировании данной математической модели.

О динамике изменения толщины и напряженного состояния оболочки можно судить по рис. 4—6, где соответственно представлены h - h(s), ain = cr^CO, a/ = аДл) --графики распределения толщины и меридиональных и окружных напряжений но длине оболочки (в безразмерном виде) для указанных выше дискретных моменте» времени. s = z/L() — безразмерная аппликата горизонтальной плоскости сечения свободной поверхности оболочки в ненагруженном состоянии, 0 < .v < 1.

Толщина оболочки h h(s, I)

Рис. 4

2006

№ 10

Окружные напряжения ос = ас0,/)

Рис. 6

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

Затраты компьютерного решения данной задачи Коши с шагом 10~4 с на интервале физического времени Т= 10 с — менее 3 минут.

В заключение отметим, что результаты расчетов для более редкой сетки дискретизации, когда п = 10, отличаются от приведенных выше для случая п = 100 менее, чем на 15 %.

СПИСОК ЛИТЕРАТУРЫ

1. Русанов П. Г. Построение расчетных моделей динамики сплошной среды с помощью метода твердых тел / М.: МАИ: сб. докл. Научн. Сим поз. «Динамические и технологические проблемы механики конструкций и сплошных сред». — 1995. —• С. 42—43.

2. Р у с а н о в Г1. Г. Анализ динамики контактного микропереключателя методом твердых тел (там же, С. 43—44.).

3. Руса н о в П. Г. Применение метода твердых тел в динамике деформируемого тела / Пермь: сб. докл. 10-той Зимней школы Ур. Отд. РАН по механике сплошных сред. — 1995. — С. 213—214.

4. Русанов П, Г. Моделирование цунами методом твердых тел (там же, С. 213)

5. Русанов П. Г., Русанова Е. М. Компьютерный анализ динамики хлопка полу кольцевой арки при симметричном нагружении силой в плоскости кольца / Подольск.: Госкомобр. 3-ий Межд. сем. «Технологические проблемы прочности». — 1996. — С. 146-148.

6. Русанов П. Г. Анализ статики пластины, нагруженной неконсервативными силами, методом твердых тел / Подольск.: Госкомобр. 8-ой Межд. сем. «Технологические проблемы прочности». — 2001. — С. 188— 197.

7. Жуков И. В., К е т а т В. В., Русанов П. Г. Анализ динамики упругой балки методом физической дискретизации // Известия вузов. Машиностроение. — 2000. —№3. — С. 3—9.

8. Русанов П. Г. Анализ динамики деформируемых тел с помощью дискретных моделей / М.: МГУ. Сб. трудов конф., посвященной 100-летию со дня рожд. чл.-корр. АН СССР Б.В. Булгакова. — 2000. — С. 116.

9. Русанов П. Г. Вычислительные модели механики на основе метода твердых тел / М.: РУ/ДН, сб. докл. на XL Всеросс. конф. по проблемам математики, информатики, физики и химии. — 2004. — С. 169—170.

10. Р у с а н о в П. Г. Примеры анализа статики и динамики пластин методом физической дискретизации / М.: РУДН, сб. докл. на XLI Всеросс. конф. по проблемам математики, информатики, физики и химии. — 2005. — С. 46—47.

П.Русанов П. Г. Компьютерный анализ нелинейных динамических процессов деформируемых механических систем по расчетным схемам на основе метода твердых тел / М.: РУДН, сб. докл. на XL11 Всеросс. конф. по проблемам математики, информатики, физики и химии. — 2006. — С. 73. 12.Руса нов П. Г. Отличия расчетов формы упругой линии гибкого стержня методами дискретизации МКЭ и МТТ // Известия вузов. Машиностроение. — 2006. — № 7. — С. 3—9.