CC BY
35
УДК 622.235.535
ОБ ОДНОМ ПОДХОДЕ К МАТЕМАТИЧЕСКОМУ МОДЕЛИРОВАНИЮ ВОЗДЕЙСТВИЯ ВЗРЫВНЫХ ВОЛН НА ПОДЗЕМНЫЙ НЕФТЕПРОВОД
А.П.ГОСПОДАРИКОВ, д-р. техн. наук, профессор, kaf_matem_spmi@mail.ru Г.А.КОЛТОН, канд. физ.-мат. наук, доцент, kaf_matem_spmi@mail.ru Е.Л.БУЛДАКОВ, аспирант, bullduckoff@mil.ru
Национальный минерально-сырьевой университет «Горный», Санкт-Петербург, Россия
Из общих уравнений механики сплошной среды, теории тонких оболочек и уравнений гидравлики составлена математическая модель нефтепровода, находящегося в скальном грунте, для дальнейшего расчета его напряженно-деформированного состояния при воздействии взрывных волн. Задача сформулирована в плоской постановке, для прямого интегрирования исходной системы уравнений выбран метод конечных разностей. На контакте массива и трубопровода рассмотрены краевые условия вида «проскальзывание» и «жесткое защемление».
Ключевые слова: математическая модель, нефтепровод, взрывная волна, скальный грунт.
В настоящее время наращивание объемов перекачки нефти обусловливает необходимость строительства дополнительных ниток нефтепроводов в одном коридоре с действующими. При прокладке нефтепроводов в скальных грунтах сооружение траншей ведется с использованием энергии взрыва, в связи с чем возникает необходимость оценки воздействия энергии взрыва на эксплуатируемый (действующий) нефтепровод.
Напряжения, возникающие в стенке нефтепровода при сейсмическом воздействии взрыва, зависят от многих причин, но, в первую очередь, от расстояния между нефтепроводом и зарядом взрывчатого вещества (ВВ), глубины заложения заряда ВВ и его массы. Для определения влияния воздействия взрыва на действующий подземный трубопровод необходимо провести большое количество натурных испытаний, что экономически и технически не всегда возможно. Поэтому для изучения воздействия взрыва на трубопровод в работе используется численное моделирование взаимодействия продольной волны, распространяющейся в упругой среде, с оболочкой, заполненной жидкой средой.
Изменение состояния упругого массива горных пород определяется системой дифференциальных уравнений в частных производных первого порядка, объединяющей в себе как уравнения движения, так и продифференцированный по времени закон Гука [1]:
дu ди
э7
до г
1 до ге 1 . - +--—+ -(огг -оее);
дг г де г
до„
1 до
- + —
ее
2
+ — о
дг г де г
ге■>
до дu (1 ди 1 , —— = — + (1 - 2Ь)|--+ — u |;
д дг I г де г
доее . дu 1 ( ди
—— = (1 - 2Ь)— + -| — + u дt дг г I де
д°-= ь[^ +1 (^-и
(1)
дt
дг г I де
Систему (1) можно записать в матричной форме:
8U 8U 1D 8U 1
-+ A-+— B-= — QU,
8t 8r r 89 r
где и = (и, и, стгг, стее, стг9)т - матрица-столбец неизвестных; Т - символ транспонирования. Постоянные матрицы А = -(ау)5х5; В = -(¿у) и 2 = -(#у) являются квадратными матрицами 5-го порядка с элементами
а13 = а25 = а31 = ¿15 = ¿24 = ¿42 = #13 = "#14 = #41 = 1 ; а41 = ¿42 = #31 = 1 " ^ а52 = ¿51 ="#52 = ¿, #25 = 2;
остальные элементы этих матриц равны нулю; и, и - радиальная и тангенциальная компо-
„ , 1 - 2у
ненты вектора скоростей; стгг, стее, стг9 - компоненты тензора напряжении; ¿ = —,-г; V -
2(1 -V)
коэффициент Пуассона.
Матричное уравнение (2) записано в безразмерном виде. Здесь координата г отнесена к характерному радиусу срединной поверхности оболочки (трубы) Я ; время I отнесено к
"ЩГУТ
величине R / c, где c =.
'р(1 + vXl - 2v)
скорость распространения продольных волн в мас-
сиве; Е - модуль Юнга; р - плотность среды. Компоненты вектора скоростей и, и отнесены к величине с, а напряжения стгг, стее, стг9 - к величине рс2.
В начальный момент времени (^ = 0) состояние массива горных пород, примыкающего к трубопроводу известно: и(0, г, е) = и0(г, е), где и0(г, е) - заданный вектор. Перед фронтом падающей волны и0 = 0; за фронтом и0 = (и0, и0, ст°г, стее, ст^У, где и0, и0 - компоненты вектора скоростей в начальный момент времени; ст°г, стее, ст°е - компоненты тензора напряжений в начальный момент времени. Значения компонент тензора напряжений за фронтом волны определяются по формулам [5]:
и0 = -ст°(^)со8 9;
и0 = a°(^)sin 9; a0r = -ст°(^)(1 - 2b sin2 9); аее = -°0(^)(1 - 2b cos2 9); a09=a0(^)b sin 29.
(3)
В формулах (3) функция ст0задает распределение нормальных напряжений за фронтом волны при t = 0 ; - расстояние от точки массива до фронта волны.
Для прямого численного интегрирования системы (1) искомую область изменения переменных (г, е), т.е. расчетную область, разобьем на прямоугольные ячейки прямыми г = гк (к = 1, 2, ..., К) и е = еу (( = 1,2,...,/), а в пространстве
(г, е, 0 выделим элементарный параллелепипед объемом V (рис.1).
Интегрируя уравнение (2) по объему V и преобразовав полученный объемный интеграл по формуле Гаусса, имеем
Рис.1. Конечно-разностная схема для определения напряженно-деформированнного состояния (НДС) массива горных пород
X
| [иС08(я,^ + лиС08(я,г) + -ви^(и,еШ5 = J-eUdV, (4)
5 ^ Г ' V
где 5 - боковая поверхность параллелепипеда; п - внешняя нормаль к ней; V _ объем параллелепипеда.
Считая, что на каждой грани элементарного параллелепипеда вектор и сохраняет постоянное значение, с точностью до малых первого порядка получим
( - ищ)е + Л (+и - й„,-) +1В (и^+1 - и^)К 1 =1б^ЛАт . (5)
В разностном уравнении (5) верхними индексами обозначено значение искомого вектора и на верхней грани параллелепипеда, крышками - его значения на боковых гранях, нижними индексами - значения на нижней грани. Для определения этих величин применим метод «расщепления» [1], в соответствии с которым значение вектора и на боковой грани параллелепипеда V находится из решения одномерной задачи Римана о распаде произвольного разрыва:
Л(ик+и - иК]) = фАиК] - ф-Аик+и - фи-1,-; вфк-+! - и) = фвик- - ФВи,,-+! - ФАи,,--!.
Здесь матрицы Ф± и ФВ определяются по формулам Ф± = Л^М±Лл ; ФВ = ЛВМ±Лв ; фл = Ф+ + Ф- ; ФВ = ФВ + ФВ; М± = -2(м| ±М), где М - диагональная матрица 5-го порядка собственных значений матриц А и В, причем матрицы А и В обладают полным набором собственных значений ( = 1,2,..., 5), одним и тем же для обеих матриц; Лл, ЛВ - матрицы собственных векторов, отвечающих собственным значениям матриц А и В.
После необходимых преобразований уравнение (5) примет вид
ик = (Е--^ФЛ --рФВ +1 +(+:,-■ + Ф+Ли,-1,-)+(фВи,,-+1 +ФВик,;-1 ).(6)
V К л гКе В г
С помощью формулы (6) осуществляется переход с временного слоя, отвечающего моменту времени t = t', на следующий временной слой t = t'+х (т - временной шаг интегрирования). Формула (6) дает возможность найти новые значения вектора и только во внутренних ячейках расчетной области, т.е. в ячейках, не имеющих общих точек с границей оболочки. Для нахождения решения (и_1,-, ик, ик, J + 1) в граничных ячейках (к = 0, - = 0 и - = J) необходимо привлечь краевые условия. Так как величины u, огг, оее четны по координате е; и; оге нечетны, то на линии симметрии е = 0,л; и = 0; оге = 0. Вводим фиктивные ячейки (к,-1) и (кД+1), полагая
и к ,-1 = diag (1,-1,1,1,-1)ик ,0;
и^J+! = diag(1,-1,1,1,-1)UkJ .
На линии контакта массива горных пород и оболочки имеем два условия, которые в
матричной форме можно записать в виде
5и-и = Ж ,
где 5 - двухстрочная матрица, зависящая от вида контакта, при «проскальзывании»
(10 0 00 ^ (10000^ 5 = , при «жестком защемлении» 5 = ; Ж - двухкомпонентный вектор; f -
V 00 - f 01 / V01000У
коэффициент трения.
Используя основные положения теории тонких оболочек, считаем, что деформации и углы поворота малы, материал изотропен и подчиняется линейному закону Гука. Принятие этих ограничений позволяет описать напряженное состояние трубопровода известными уравнениями Новожилова в рамках плоской деформации и с учетом продифференцированного закона Гука [2]:
, 8и 8T
ph ä= ш +N + 90;
, 8w 8N ph— =--T + qr;
8t 89 r
, 3 8ra 8M ph3
8t 89 8T Eh (8u 8t " 1 -v21 89
- N;
+ w l;
(7)
8M
Eh3
8ra
дt 12(1 -V2) де
дN , <: дw -= пи\ юн---и I.
дt у де
Система дифференциальных уравнений первого порядка в частных производных (7) в матричной форме имеет вид
8V 8V
— + F— = CV + P, 8t 89
(8)
где V = (и,w,ю,Т,М,- вектор неизвестных; ^ = -(Уу)6х6,С = (сгу)6 ,Р = (ргу)6х1 - постоянные матрицы 6-го порядка с элементами
У = г = У = =" =" = г = г = = Еп =" = ; с = ;
У14 = У 26 = У 35 = С1б = С24 = С3б = , ; У41 = У53 = С42 = , 2 ; у62 = С61 = ; С63 = , ;
рп 1 - V п
#е #г
р- =-т; Рп = -т, рп рп
остальные элементы этих матриц равны нулю; и и w - радиальная и тангенциальная компоненты
„ _ п2 t'
вектора скоростей; ю = — ю ; ю - скорость сдвига; И - толщина оболочки; Т и N - касательные и г+х нормальные усилия соответственно; М - изги-
Е г бающий момент; G =-; Е - модуль Юн-
2(1 + V)
га; V - коэффициент Пуассона; р - плотность материала трубы; #е и #г - компоненты поверхностной нагрузки.
Уравнение (8) записано в безразмерном виде. Здесь плотность материала трубы р отнесена к р0 - плотности материала массива; компоненты вектора скоростей и, w отнесены к скорости продольной волны в массиве С0; усилия N и Т - к величине Яр0С^; время t -
40
Vi
F Q V+1
Vi
i * h9 * i+1 9
Рис.2. Конечно-разностная схема для определения НДС оболочки
к Я/с0; изгибающий момент М - к Я2р0с0; скорость сдвига ю отнесена к с0 / Я; толщина трубы п отнесена к ее радиусу Я; параметры Е ,0, #е, #г отнесены к р0 с0.
Интегрируя матричное уравнение (8) по объему (прямоугольник) О = [^ t + т]х[ег, е^] (рис.2), получим
| (FVdt - Vd9) = Л (^ + Р)& ёе ; (9)
V = V, - -1 F(Пг- П) + А т
пе
(10)
где Вг = (СV + Р); П,+1 = ^; П, = .
t t
Здесь верхним индексом обозначено значение вектора неизвестных V на верхней стороне прямоугольника, крышкой - его значения на боковых сторонах, нижним индексом -значения на нижней стороне. Для определения этих величин, как и при расчете вектора и для массива, применяется метод «расщепления» [1]:
т т
V = (Е - — Ф)Г, + — (Ф+V- + Ф-ГМ) + т
п9 п9
т т
(Е-—Ф) Д + — (Ф+Д-1 + Ф - им) е 2пе
(11)
Для определения векторов V--, 0-1, V/+-, В1+1 используем симметрию относительно лучей е = 0, е = % , полагая в фиктивных ячейках с номерами -1 и I +1 следующие значения:
V-- = diag(-1,1,-1,1,1,-1)^ ; VI+1 = diag(-1,1,-1,1,1,-1)^ ; Р-1 = diag(-1,1,-1,1,1,-1)Р0; Р+1 = diag(-1,1,-1,1,1,-1)Р1; Б-1 = СУ_Х + Р-1; Д+1 = ^+1 + Рг+1 .
Движение жидкости внутри оболочки описывается известными линеаризованными уравнениями Эйлера и уравнением неразрывности с учетом сжимаемости жидкости [3]:
ди„
дt
дие
дt
др; дг' 1 др ; 7 дё'
(12)
др 1 д . ч 1 дие
— =---(игг)---е.
дt г дг г де
Система дифференциальных уравнений (12) в матричной форме имеет вид
дА тдА дА — + Ь— + С— = КА дt дг де
(13)
где А = (иг, ие, р)Т - вектор неизвестных; Ь = (¡у )3х3, С = С ) К = (ку )33 - матрицы 3-го порядка с элементами ¡13 = ¡31 = 1 ; с23 = с32 =-к31 = 1/г, остальные элементы этих матриц равны нулю; г, е - текущие координаты; р - давление; иг, ие - компоненты вектора скоростей.
Матричное уравнение (13) записано в безразмерном виде. Здесь компоненты вектора скоростей иг, ие отнесены к скорости распространения продольной волны в жидкости а0; время t отнесено к Я/а0; давление р - к р^ . Решение уравнения (13) осуществляется тем же методом «расщепления» [1], использованным при решении уравнения (3), основное отличие заключается только в наборе собственных значений матриц Ь и С.
Расчетная формула, описывающая изменение состояния жидкости, имеет вид
г
AkJ =
тт E--ФL--Фс + Кт
V hr he J
AkJ +(ф LA k+1,J + Ф+Л-1,J )+ (ФCAk,J+1 + Ф C Ak, J - ). (14)
Формула (14) дает возможность найти новые значения вектора А только во внутренних ячейках расчетной области, т.е. в ячейках, не имеющих общих точек с границей оболочки. Поэтому необходимо ввести краевые условия на линии контакта оболочки и жидкости, которые могут быть записаны в матричной форме:
cSA-1,J = WaO, (15)
где S - однострочная матрица, имеющая вид S = (l, О, О ).
Таким образом, замкнутая система дифференциальных уравнений (6), (11) и (14) с учетом граничных и начальных условий представляет собой математическую модель совместного движения грунта, трубопровода и заполняющей его жидкости в условиях внешнего динамического нагружения. Сформулированная краевая задача на основе разработанного численного алгоритма сведена к вычислительной программе на алгоритмическом языке Фортран-9О [4], с помощью которой можно оценить влияние воздействия сейсмовзрывной волны на напряженно-деформированное состояние трубопровода.
ЛИТЕРАТУРА
1. Валландер С.В. Лекции по гидроаэромеханике. Л., 1978. 296 с.
2. Годунов С.К. Разностные схемы / С.К.Годунов, В.С.Рябенький. М., 1973. 4ОО с.
3. Новожилов В.В. Линейная теория тонких оболочек / В.В.Новожилов, К.Ф.Черных, Е.И.Михайловский. Л., 1991. 656 с.
4. Brian D. Hahn, Edward Arnold. Fortran 9O for Scientists and Engineers. University of Cape Town, 1994. 198 p.
5. Nelson I., Baron M.L., Salander I. Mathematical models for geologic materials for wave propagation studies // Shock waves and mechanical properties solids. Syracuse: Univ. press, 1971. P.289-351.
REFERENCES
1. Wallander S. V. Lekcii po gidroaeromekhanike (Lectures on hydromechanics). Leningrad, 1978, p.296.
2. Godunov S.K., Ryaben'kii V.S. Raznostnie skhemi (Difference schemes). Moscow, 1977, p.400.
3. NovozhilovV.V., Chernykh K.F., Mikhailovsky E.I. Lineynaya teoriya tonkikh obolochek (Linear theory of thin shells). Leningrad, 1991, p.656.
4. Brian D. Hahn, Edward Arnold. Fortran 90 for Scientists and Engineers. University of Cape Town, 1994, p.198.
5. Nelson I, Baron M.L., Salander I. Mathematical models for geologic materials for wave propagation studies. Shock waves and mechanical properties solids. Syracuse: Univ. press, 1971, p.289-351.
MATHEMATICAL MODELING OF THE IMPACT OF BLAST WAVES ON UNDERGROUND PIPELINES
A.P.GOSPODARIKOV, Dr. of Engineering Sciences, Professor, kaf_matem_spmi@mail.ru G.A.KOLTON, PhD in Physics and Mathematics, Associate Professor, kaf_matem_spmi@mail.ru E.L.BULDAKOV, Postgraduate student, bullduckoff@mail.ru National Mineral Resources University (Mining University), St Petersburg, Russia
Mathematical modeling of the impact of blast waves on the underground pipeline was composed from general equations of continuum mechanics, shell theory and hydraulics equations. The problem is formulated in the plane formulation for direct integration of the native system of equations' chosen method of finite differences. At the contact of the array and the pipeline, boundary conditions of slippage and rigid clamping are considered.
Key words: mathematical model, oil pipeline, blast wave, rocky ground.
