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

№ 3 (45) 2013

И. В. Бугай, аспирант Московского государственного технического университета им. Н. Э. Баумана

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

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

введение

Численное моделирование последствий воздействия потоков излучений различной физической природы на тонкостенные несущие элементы конструкций включает в себя теплофизические, термомеханические, газодинамические и прочностные аспекты [1-4].

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

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

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

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

Уравнение переноса энергии твердой фазы с учетом ФХП и переноса тепла газовым потоком продуктов разложения связующего записывается для каждого подслоя материала. Система этих уравнений с соответствующими граничными и начальными условиями решается методом конечных разностей по неявной схеме [5]. Полученная система нелинейных алгебраических уравнений решается итерационным методом Ньютона, а на каждом шаге итерации для системы линейных уравнений с трехдиагональной матрицей применяется метод прогонки [5]. Для контроля устойчивости и точности счета в процессе вычислений проверялось выполнение закона сохранения энергии. Для проверки правильности численного алгоритма и программ, реализующих модель нагрева и теплового разрушения композитных материалов при многократном воздействии потоков излу-

№ 3 (45) 2013

1

и

э

I &

1

0

1

!

I е

о §

I

I

§

со §

о

§

£5 *

ч

I

со

0

1

I

чении, проводилось тестирование на ряде задач, имеющих аналитическое решение [6].

Механическое действие. Для расчета параметров действия потоков излучении допустимо использование одномерной плоской модели [3], которая основана на численном решении нестационарных уравнений механики сплошной среды. Используемая система уравнений движения, траектории частиц, неразрывности и энергии позволяет единым образом описать процессы в конденсированном и газообразном состояниях вещества. Эта система уравнений дополняется широкодиапазонными уравнениями состояния [7], определяющими уравнениями деформирования и разрушения, а также моделями пористости [8]. Особое внимание при реализации программы расчета волновых процессов было уделено ее адаптации к определяющим уравнения состояния гетерогенных защитных покрытий, построенным методом элементарной ячейки [9].

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

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

Уравнения движения оболочки вращения записываются в виде (1) (в — координата вдоль образующей поверхности приведения; ф — угловая координата) [11].

Напряжения в слое ортотропной оболочки на расстоянии г от поверхности приведения записываются в виде (2) [12].

дМ)- СЯ^+кт=^

Эв Эф Св ф 1 8 вЬ дt2

ач+а^+^ в+^ ЯОф=т^ *

Эф Эв Э(ЯОв) ЭОф

- я (+к ыф)+р(в,ф,о=

Эв Эф

=твЯ-:-г,

_ = 2 д(яМв)+2М- СЯ.м

в Я Эв Я Эф ЯСв ф,

о 1 Э(ямвф) + М

Ч =Я Эф + Я дв Яс1вМфв,

, 1 С9 . 1 ^(9)

к,=—=— к2 =—=-—,

1 Я Св 2 Я Я

(1)

9=аг^д

(С2 X / Св2 ^

С2 Я / Св2

где Nф, Ывф, N — усилия, М„ Мф, Мз¥ Мфв — моменты, Ов, Оф — перерезывающие силы в поверхности приведения, отнесенные к единице длины дуг соответствующих координатных линий; Я(в) — радиус кривизны параллельного круга; 8 — угол между радиусами кривизны нормального сечения и параллельного круга; X = Х(в), Я = Я(в) — параметрическое задание образующей оболочки вращения, р(в,ф,^ — давление механического действия потоков излучения.

о в =

Е

Оф=-

1—V2—д2

Е

1—V2—д2

о вф=оф в =

[(1+д)(е в+к вг)+\<Еф+Кф z)], ■[(1—д)(еф+кфг)+^Ев+к(2) [со+х г],

2(1+у)

где приведенные модуль Юнга Е, коэффициент Пуассона V и поправочные коэффициенты ■&, ^ определяются по известным характеристикам ортотропного материала Е1, Е2, П12, П21, ^12'

92

№ 3 (45) 2013

^2 V21 V

Е, "Е2 Е'

1 1-;

Е = Е ,

1= 1+;

Е2 = Е ,

1 2(1+¥)

а12 =Е

откуда

Е = 2 Е E2

V =

Д = -

E1 + E2'

E2 V12 + E1V21 Е

E1 - Eг Е + E2,

(3)

¥ =

2а,

—1.

Из (2) следует, что параметры V, Д, ¥ должны удовлетворять неравенствам:

V2 + Д2 < 1, ¥ > -1.

Используя (2) и интегрируя по толщине многослойной оболочки, находим:

N = [о 8 (1 + к2 z)dz = = в+е 8 + Д* еф + 5Х + кф, ^ = [оф (1 + к^ =

= + Вд£ф+ + 5ДКф ,

М3 = [о 8 (1 + к2 z) zdz =

= + ^2еф + + ^2Кф , Мф = |оф (1 + k-z) zdz = = $ ие0 + 5-еф + Оик + D;кф,

¥2 т,

^* = [о 8ф (1 + кг z)dz = В¥2 ® + N 8 = [о 8ф (1 + k1z)dz = В¥1Ш +

Мэф = [о8ф (1 + к2z)^ = 5¥2® + ^2т, Мф 8 = [о 8ф (1 + k1z)zdz = 5¥1Ю + D¥1т,

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

Отметим, что усилия и моменты из (5) с жесткостями (7) тождественно удовлетворяют уравнению равновесия моментов относительно оси z:

Чф- N 8 + кМф- к2 Мф 8 = 0.

Входящие в соотношения (4), (5) параметры деформаций определяются по полям перемещений [11]:

ди , г8 = - + ^

£ф =

дv dR

- + и + кгШ,

Ядф Rds

(4)

(5)

ди п д ( V

® =-+ Я— I —

Ядф дs I, Я

д2w dk1 , 2

К = -—т + —1 и -s дs

Кф = --

ds дw

Ядф ^ Ядф

д

dk1 2

+ —1 и - k;w, ds 2

(6)

2

т=-

г

( д2w dЯ дwЛ

дsдф Rds дф

+(к - кг)

1 ди - Яд( V

Я дф дs I, Я

Начальные условия для системы уравнений (1), (4-7) записываются в виде (8).

£

из со 55

Вд+ = [ 1Е(12+Д^2 (1 + к2z)dz, вД = [ 1Е(12+ДД2 (1 + к2z)zdz, DД = [-

1 - V2 - ;2 ■' 1 - V2 - ;2 ■' 1

Е(1 + Д)

-4—М1 + к2z)z dz,

-V2-;

- V2 - Д2(1 + = [+ zdz, ^ = [^2 - Д2 ......

Е(1-Д)

в; = [ Т

= [ (1 + к1^, ^ = [ (1 + ^zdz, ^ = [ 1

Е (1-Д)

(1 + k1z)z гdz,

В^2 = [ 1 Е2У а2 (1+ к2z)dz, ^ = [ 1 ^ ^ (1 + ^Z)ZdZ, Dyг = [" 1 - V2 - ;2 1 - V2 - ;2 ■' 1

Е^(1 + k1z) z

- V2 - ;2

Е ^1+;)

(7)

2 2 (1 + к2z)z dz.

- V2 - ;2

93

№ 3 (45) 2013

du dt

и\

dv dt

= v\

= w

= 0,

= 0, ^ dt

Ip (s)

= V0(s) = , (8)

t=0 msh

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

В качестве граничных условий при ф = ±п /2 требуется выполнение соотношений симметрии:

= 0, —

=±п/2 Эф

=±п/2

= 0, ^ Эф

= 0,

=±п/2

d3w Эф3

(9)

= 0.

±п/2

1

I &

» I

о

s

! Ü

0

1

и

I

§

со

s

о

§

ÎS *

4

I

s 00 О

!

5

На краях оболочки задаются следующие основные граничные условия [11]: а) свободный край

Ns 0 = 0, Ms 0 = 0,

sls=0 sls=0

Ns ф+ k2Msrn\s=_0 = 0,

(10)

а + 1 dMs ф

R Эф

= 0;

=±п/2

б) шарнирно-закрепленный край

M\ п = 0, и 0 = v 0 = w 0 = 0; (11)

sis=0 ' ls=0 Is=0 ls=0 ' 4 '

в) абсолютно-заделанный край

dw .

--+ ки

ds 1

= 0,

(12)

и n = v n = w п = 0.

Is=0 ls=0 ls=0

конечно-разностные соотношения для расчета деформирования ортотропной оболочки вращения переменной толщины

Явная конечно-разностная схема строится аналогично одномерной схеме, рассмотренной в [13]. Перерезывающие силы Qs, Оф и усилия определяются в серединах

соответствующих сторон, перемещения и моменты — в узлах, моменты Мвф, Мфв — в центрах ячеек (рис. 1). Уравнения движения в конечных разностях имеют следующий вид:

d2U ^ 2(Ri+V2Nsi+112 j -Ri-M2Ni-M2 j )

j

'/+1/2 s/+1/2j ''/-1/2' "/-1/2j /

(As +As,J R

2(N

-Nn )

фs//+V2 /чфs/j-1/2 / , T^n + k1/Qs/j

(Аф , +ф j-1) R/

dR ^ Nn.,

ф j

Rds

m

2n

d2v

sh 2

dt2

j

_2(n;7+112 - n;i1 -1/2 )

" (Аф j +ф H) R l

+ 2(R/+1/2Nsф/+1/2^ R/-1/2Nsф/-1/2j )

(As, tAs,.,) R, k2 ] N«,

/

d2w Y 2(R(+1/2 Q;

^Tsh "dt2 y .

V

_R Qn )

7+1/2 ^s /+1/2 , /-1/2 s/-1/2 , J

2(аф// +1/2 Q<n/j-1/2 )

(Аф j +ф j-/) Rl

(As+As/-1) R

- k1 Nsn - k2/ N,+p

(13)

Qn = R MMnsi+1j -ям:,

^ = R/+1/2 As/ +

фs/+1/2 j+1/2 -Mi>s/+1/2 j-1/2 )

2(M

(Aф j+Aф j-1) R,

! Mn Rds I ,

Q

m;,+1 -mi,

ф/, +1/2 "

R Aф ,

2(Msф(+1/2 j+1/2 Msф /—1/2 ,+1/2 )

(As+As,-1) R

dR ) Mn

Rds I. MфS(j+1/2,

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

94

t=0

v

s=0

№ 3 (45) 2013

Qsij = 1 (/+1/2 I + -1/2 I ), Nsj =1 ((+!/2 I + _п I),

Чф I=1 ((

МП+1/21 = 1 ((

+ Nn

sф / +1/2 I т '^ф /-1/2 I

'/+11 + Кп)

Оф ,= 1 ( NФ ! = ^ N

+ Оп )

■ /1+1/2 т ^ф //-1/2 )'

N=1 (

+ Nn )

/1 + 1/2 чф //-1/2 ,/'

Мф

s/j+1/2

2

= 1 (

+ Nn 1

s/1+1/2 т ' чф s/1-1/2 ,)

+ Мп )

в/+11+1/2 т 1 "'ф s/1+1/2

(14)

£ иа

со

55

Параметры деформации (7) определяются в тех узлах, в которых они требуются для определения усилий или моментов по соотношениям (4), (5), в соответствующих точках определяются и жесткостные характеристики (6). Например, для моментов М8ф, Мф8, определенных в центрах ячеек, получаем:

пп

М 1+112,+112фв = °¥1/+1/21+1/2 ® / +1/21+1/2 +

+D тп

т|^¥1/ +1/21+1/2 1 /+1/21+1/2 '

ю

и" + и - и - и

и/+11+1 ^ uij+1 и!+1/ и/7

/+1/21+1/2

Я

/+1/2

Я/ +1/2

V

/+1,1+1

2Дф(

2Дs,

Я+

V". 1 V"1• ', I+1 + 1+1,1

ЯЯ+

Я

(15)

(16)

и, У,

' У

э' Эф' э

Рис. 1. Расчетная сетка

т

.-wn - wn < + wn

/+1/+1 /+11 //+1 /1

i+1/2j+1/2 "

Я

/+112 ,п

/+7+1

V

Дф,

' dR ^ wn+1- wL, + wIní+1 - wn Л

4+11

' /у+1

, Яds

/+1/2

2Дф I

(17)

( /+112 " " к2/+1/2 )

I и/+7+1^и /7+1 и/+7 ип

Я

/ +1/2

2Дф,

Я

/+1/2

V

/+1,/+1

2Дs,

Я+

V" 1 V'11 • /,/+1 + /+1,/

Я Я+ Я

' У

ДГ < min

ДГ < min

Я Дs Дф

^ дs ^ЯДф

В; \ В- ^

в =

Дs

Я Дф'

(18)

Рассматриваемая конечно-разностная схема имеет второй порядок аппроксимации по времени и пространственным переменным на равномерных сетках. Используя критерий Неймана [10], для рассматриваемой конечно-разностной схемы получаем практические условия устойчивости в виде (здесь минимум берется по всем расчетным ячейкам) (18):

Результаты расчетов совместного теплового действия излучений и нестационарного давления

Для демонстрации работоспособности предлагаемого комплекса программ далее представлены результаты расчетов совместного действия физических факторов различной природы (солнечной радиации (СР), объемно-поглощаемого излучения (ОПИ), светового излучения (СИ) и ударной волны (УВ) мощного взрыва [14]) на ортотропную свободно-опертую по контуру панель малой кривизны. Геометрические размеры пане-

95

1

п

п

+

2

п

№ 3 (45) 2013

— т,°с д = 0,6Вт/см2 + СР

— Я = 0,1Вт/см2 + СР

д = 0,025Вт/см2 + СР

О 0,2 0,4 0,6 0,8 1

Рис. 2. Стационарное распределение температуры при совместном действии ОПИ и СР

ли: Я1 = Я2 = ^; h = 4 мм, Lx = 0,5 м (длина); § Ly = 0,17 м. Панель изготовлена из высоко-5 температурного стеклопластика, допускаю-|| щего нагрев до Т = 250°С без существенной ® потери жесткости и прочности. £ Нагрев стенки постоянно действующими факторами (СР и ОПИ) ограничивается кон-

5 векцией и тепловым излучением внешней Ё, поверхности стенки. Коэффициент конвекции принимается равным а = 10-3 Вт / см2/К"

Е, [15]. Результаты расчетов устанавливающегося при воздействии стационарно-

§ го профиля температур для потоков ОПИ

| с q = 0,1 Вт/см2, 0,6 Вт/см2, 0,025 Вт/см2 и СР

® с q = 0,1125 Вт/см2 представлены на рис. 2.

| При воздействии одной лишь солнечной ра-

| диации температура постоянна по толщине

¡ч стенки и может достигать величины Т= 70°С.

* Совместное действие СР с ОПИ приводит

| к практически равномерному по толщине

§ профилю. В этом случае нагрев больше,

§ и при q = 0,6 Вт/см2 температура может дос-

6 тигнуть величины Т= 135°С, что не приводит к существенному изменению деформацион-

| ных и прочностных свойств рассматривае-§ мой панели.

о

^ Еще больший нестационарный нагрев конструкции реализуется при воздействии ?! светового излучения мощного взрыва. Ре-

зультаты расчета нагрева стенки панели световым излучением с учетом (жирные линии) и без учета постоянных сопутствующих факторов (СР и ОПИ с q = 0,6 Вт/см2) представлены на рис. 3, 4 для W = 440 Дж/см2, ттах = 1 с и W = 15 Дж/см2, ттах = 0,3 с для моментов времени, соответствующих окончанию облучения СИ и приходу УВ ^ = 8 с).

О б л уч ен и е С И с п арам етрам и W = 440 Дж / см2, ттах = 1 с приводит на момент прихода УВ к нагреву стенки по всей толщине до Т = 300-400°С, что недопустимо для стеклопластиковых материалов. В связи с этим исследования действия УВ с соответствующими параметрами (Арш = 0,75 атм, Dsw = 435 м/с) не проводились.

Облучение СИ с W = 15 Дж/см2, ттах = 0,3 с приводит на момент прихода УВ к прогреву стенки по толщине до Т = 160-170°С, и для этого варианта воздействия имеет смысл исследовать панель на прочность к действию ударной волны (Ар^ = 0,09 атм, Dsw = 350 м/с). Считаем, что УВ набегает вдоль более длинной стороны панели, а нормаль к поверхности УВ параллельна поверхности панели. При этом давление на стенку

Рис. 3. Нагрев СИ ^ = 440 Дж/см2)

№ 3 (45) 2013

350 - т,° с

300 -

250 -

200 - 1 = 8с

150 - 1 = 0,6с

100 - \. = 8с

50 -

О 0,2 0,4 0,6 0,8 1

Рис. 4. Нагрев СИ ^ = 15Дж/см2)

меняется по закону (0(х) — тета-функция): р^, х, у) = Држ - х).

Модули упругости и коэффициенты Пуассона при Т = 20°С принимаются равными Е1 = 27 ГПа, Е2 = 24,3 ГПа, у12 = 0,3, у21 = 0,27.

Рис. 5. Прогиб по х (у^у/2)

На рисунках 5, 6 представлены результаты расчетов профилей прогиба м для моментов времени t = 0,1; 0,2; 0,3; 0,4; 0,5; 1; 2; 55 3; 4; 5 мс. Видно, что прогиб не превышает 1/5 толщины оболочки, что вполне допустимо для условий функционирования рассматриваемой панели.

На рисунке 7 показаны изменения во времени максимальных (по всей панели) деформаций (тонкими линиями показаны максимальные деформации панели без учета нагрева). Поскольку величина этих деформаций не превышает 0,06%, то совершенно очевидно, что рассматриваемое совместное тепловое действие потоков различной физической природы (СР, ОПИ и СИ) и механическое действие УВ волны не представляет опасности для исследуемой панели.

Заключение

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

Рис. 6. Прогиб по у (х=^/2)

№ 3 (45) 2013

0,06

0,04

0,02

l«Wl. %

1

I

с® &

I

0

1

t

1

б

0

1

n

I

oo

Si

о

§

i»

t

4

I

s

CO О

!

5

0 1 2 3 4 5

Рис. 7. Максимальные деформации

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

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

Список литературы

1. Острик А. В., Слободчиков С. С. Математическая модель разрушения композитных оболочек высокого давления под действием лучистых потоков энергии // Математическое моделирование. 1995. Т. 7. № 10. С. 33-46.

2. Горшков А. Г., Дергачев А. А. Воздействие высокоинтенсивного потока энергии на элементы конструкций из композиционного материала // Механика композитных материалов и конструкций. 1996. Т. 2. № 1. С. 51-68.

3. Грибанов В. М, Острик А. В., Слободчиков С. С. Тепловое и механическое действие рентгеновского излучения на материалы и преграды // Физика ядерного взрыва: в 2 т. Т. 2. Действие взрыва. М.: МО РФ ЦФТИ, 1997. — 256 с. С. 131-194.

4. Грибанов В. М, Острик А. В., Ромадинова Е. А. Численный код для расчета многократного комплексного действия излучений и частиц на многослойный многофункциональный гетерогенный плоский пакет. Черноголовка: ИПХМ РАН, 2006. — 92 с.

5. Годунов С. К., Рябенький В. С. Разностные схемы. Введение в теорию. 2-е изд. М.: Наука, 1977. — 440 с.

6. Лыков А. В. Теория теплопроводности. М.: Машиностроение, 1984. — 264 с.

7. Бушман А. В., Ломоносов И. В., Фортов В. Е. Модели широкодиапазонных уравнений состояния при высоких плотностях энергии. Препринт МВТ АН СССР. М., 1990. № 6. С. 287.

8. Херрманн В. Определяющие уравнения уплотняющихся пористых материалов // Механика. М., 1976. № 7. С. 178-216.

9. Острик А. В., Острик Е. А. Квазистатическая модель установления давления в многокомпонентном пористом гетерогенном материале при воздействии излучения // Математическое моделирование. 2001. Т. 13. № 6. С. 93-98.

10. Рихтмайер Р., Мортон К. Разностные методы решения краевых задач. М.: Мир, 1972. — 418 с.

11. Амбарцумян С. А. Общая теория анизотропных оболочек. М.: Наука, 1974. — 448 с.

12. Работнов Ю. Н. Механика деформированного твердого тела. М.: Наука, 1988. — 712 с.

13. Witmer E. A. [et. all] Large dynamic deformations of beams, rings, plates and shells // AIAA J. 1963. № 1 (8). Р. 1848-1857.

14. Верховенцев В. Н., Лоборев В. М., Фортов В. Е. [и др.] Физика ядерного взрыва: в 5 т. Т. 1. Развитие взрыва. М.: Физматлит, 2010. — 832 с.

15. Основы теплопередачи в авиационной и ракетно-космической технике / под общ. ред. В. С. Авдуевского, В. К. Кошкина. 2-е изд., пе-рераб. и доп. М.: Машиностроение, 1992. — 528 с.

98