Асимптотическое представление и численный расчет конечных деформаций криволинейного подземного трубопровода Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 11, № 1, 2006

АСИМПТОТИЧЕСКОЕ ПРЕДСТАВЛЕНИЕ И ЧИСЛЕННЫЙ РАСЧЕТ КОНЕЧНЫХ ДЕФОРМАЦИЙ КРИВОЛИНЕЙНОГО ПОДЗЕМНОГО ТРУБОПРОВОДА*

О. П. ТКАЧЕНКО Вычислительный центр ДВО РАН, Хабаровск, Россия e-mail: olegt@as.khb.ru

A geometrically nonlinear 2D mathematical model of the pipeline as shells in the viscous medium is constructed. It is shown that without loss of generality the governing system of equations can be reduced to 1D problem by a special asymptotical representation of the solution. Software package for the numerical analysis of the model is implemented and the numerical calculations have been carried out.

Введение

Трубопровод, нагруженный потоком жидкости и помещенный во внешнюю среду, является сложной механической системой, поведение которой исследовано не полностью. Обзор работ по этой проблематике является темой отдельной статьи. Мы приведем примеры работ, относящихся к нашей теме и описывающих нелинейные системы. В классических работах [1, 2] заложены основы исследования нестационарного течения жидкости в трубопроводах. В монографии [3] исследовалось деформирование оболочек в потоке жидкости и газа. Статьи [4-6] и многие подобные работы посвящены изучению потока жидкости в трубке, моделирующей кровеносный сосуд. Работа [7] отличается тем, что в ней рассматривается совместное движение потока жидкости и трубы (в стержневом приближении) в приложении к трубопроводным системам.

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

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований и Администрации Хабаровского края (грант № 04-01-97004).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.

1. Физическая постановка задачи

Пусть трубопровод проложен в сильно вязкой среде и имеет форму слабо изогнутой плоской кривой Г = {х0,у0 ■ х0 = х0(в), у0 = у0(в)}. Он заполняется стационарным потоком жидкости, который движется под действием постоянного перепада давления. После этого трубопровод начинает двигаться, так как теряет равновесие.

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

2. Формулировка двумерной математической модели

Введем глобальную декартову систему отсчета {Охуг} и криволинейную лагранжеву систему координат {Ов9Я} (рис. 1) [9]. Пусть ш3, Шв, 'и — компоненты вектора перемещения стенки трубы как трехмерного упругого тела. В [8] показано, что в данной задаче можно считать конечными только деформации, отражающие связь прогиба осевой линии и ее растяжения. Тогда интересующие нас компоненты тензора деформаций имеют вид

1

\Z~9ss

дш3 дв

+ швк сое 9 + шик вт 9 —

1

2v/gSS

дшв \2

дв

+

дш

и

дв

21

дшв + щи Яд9 + Я'

дш

£ии

и

£вв

К-

1

дЯ 2 \ дв

Здесь д33 — компонента метрического тензора.

дшв + дш.

шя

Яд9 у/дЩ

к сое 9

(1)

£

Рис. 1. Система координат.

Пусть и, V, т — компоненты вектора перемещения срединной поверхности стенки трубопровода в актуальной конфигурации (по терминологии [10]). Для перехода к оболочке воспользуемся соотношениями [11]:

т3

тв тп

/-, , ч 7 дт

и(1 + ^17) - Ая-' А д-

/-< 7 4 7 дт -(1 + к27) - в

(2)

т,

7=я - д0, в=д0, к2 = = р +8"пи

Я0 р0 + Я0 в1п и

А = 1 + к0Я0 в1п и,

где Я0 — радиус трубы; к0 = к0(з) — начальная кривизна осевой линии Г; р0(з) — начальный радиус кривизны Г.

Следуя методике, изложенной в [11], в работе [8] был проведен переход от уравнений движения трубопровода как трехмерного упругого тела к уравнениям для оболочки. После перехода к перемещениям срединной поверхности стенки получена система уравнений

2д2и' 1 — ^д2и' р£Я2^2 д2и' 1 + V д2^' дт' . ,/1 — V .

а + — "д^ — "д^ + — + ^ + е/81ПЧ~и —

+

2 д2 и' дт'

—2а ^^^ + а(1 — V)

1 — V ди' \ 3.

— а3(1 — 3е/ в1п(

1 + V д V \

■а

2 -дсд^; +е/ е°8 ч~+

3 — V дг>'

ди

1 — V 2 дV д2^' д2^' 1

-а

+

~дт' д2т' дг>' д2V'"

+ 5С2"-2и*^ сов и

-Ф^о);

д(2 ди2 Е * дт2 Е *к*

ЯЛ О, 5 — 1п

1Рдт и

Я0

Е *к

1 + V д2и +---— а

+

дт'

дт'

+е/ ет и ( — — (1 — V)а2~ос2

д2^' 1 + V д2и' \

2

а

2 д(ди ди

дг>' 3 — V ди'

дсд^; + с°8 ^+ж ——^)—

—а2(1 — 2/ в1п(

дт' д2т' дг>' д2,у' \ 2

"д^дсдй + "дсдс^ + ае/ с°8 и

дт' \2 / д-у'4 2

р4Я2д2т'

к*2 / 0д2т' д2т'

+

Е* д 4т'

+ т' + 12 (^"д^ ди' дг>'

+

ди2

+ а4

, д4т'

дС

+ 2а2

д4т'

д( 4 д(2ди2

ди' дг>' \

9 4т'\ ди' дг>' / ' ди' дг>'\

+^ + ж+е/ 81п Ч +(1—V )а^с + "ы)+

а

с°в и ■ V ——(1 — е/ в1п(

дт'

5С

+ ( ас)

Е * к

"(Р — Ре ).

(3)

Здесь £ = ^/1; т = и' = и/Я0; V = ^/Я0; т' = т/Я0 — безразмерные длина дуги, время и искомые функции; I, ^ — характерные масштабы координаты и времени; е = Я0 ■

к(с Т)

шах|к0(£)| << 1 — малый параметр; /(£,т) = - ' — функция текущей кривизны

шах|к0(( )1

оси; а = Я0/1 — коэффициент; «(£, т) — текущая кривизна осевой линии; р — внутреннее

1

О

+

2

1

давление; ре — давление внешней среды; Фг — сила внутреннего трения; и* поперечного движения трубы.

Уравнения (3) необходимо дополнить формулами [8]

скорость

Ре

РдгЯК{1 - СОв в) +

2^и* вт в

Яо 0,5-1п

7 Рдг и

4 ^

По

Р = Ра + - С) + еМУ1о ^ в,

а

Фг^о) = 2 ^вУ

2

во"

(4)

Здесь ув0, Рf — скорость и плотность внутреннего потока жидкости; ра — атмосферное давление; рдг — вязкость и плотность внешней среды; 7 — число Маскерони; По — глубина закладки трубопровода; в — коэффициент внутреннего трения [12]; С — длина трубопровода в единицах I.

Модель (3), (4) дополним однородными начальными и краевыми условиями.

3. Приведение математической модели к одномерному виду

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

Разложим искомые функции в ряды по малому параметру е:

и' = и0 + еи1 + 0(е2), у' = у0 + еу1 + 0(е2),

и;' = ио + еи>1 + 0(е2), и* = еи1 + 0(е2), (5)

и воспользуемся обычной техникой асимптотического анализа, составляя уравнения для всех порядков по е. Для того чтобы считать у0, и0 не зависящими от в, необходимо положить Е0/к0 ^ 1, тогда у0 = 0, и0 = и0((), и0 = и0((). Это означает, что мы пренебрегаем высокочастотными колебательными процессами.

Уравнения нулевого приближения по е принимают вид

а

д 2 и0

+ иа

ди0 3 / ди0 д2и, — а

д(

П*2

ио + 12 1а д(2 1

2д 2и0 , 4

+ а

д( д(2 д4и 0 ди0

1а

Е*П* 2

£ву1,

д(4

а ди 0

+ "а1к - т( ж

Е*П*

Ра + 1вУ2зо(С - I) - РдгдПо

(6)

Выражение для скорости поперечного перемещения осевой линии получим, дифференцируя по времени перемещение этой линии в направлении нормали к ней и учитывая, что У0, и0 не зависят от в:

2п

* еКош [ (ду1 ди1 . . *

и =- / -7— соя в + —— вш в ав, щ

2п ] \дт дт ' 1

0

2п

со8 в+^ - в) ав. (7)

2

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

« = Й1((,т)8т 0, VI = ^1((,г)ео8 0, и = й^^т^т (8)

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

,д 2и1 1 — v_ 1 + V

2- -1 - - ^ | - —1 . ди1 . , — —« — ~ аЦ£ + ^ + 7

1 - V

-и0 — 2а

,д2м0

5С2

+ а(1 — V)

ди0

—а

ди 52и0 + ди0 д2и

д( д(2 д( д(2

, „ з, ди 52и р4Я^2 52«1 + 3а /- —

д( д(2

1 — V 2 д^ _ 1

-а ^— v1 —

дС

2

Е *к*

Д0 0.5 - 1п

1Рдт еЩ

Е* дт2 '

1 + V ди1 + —— а-^" +

дС

, - 3 — V д«0 \ 2 +и1 + / ( и--— а— I — а

ди0 ди1 р4Я^2 д^ ^^^^ = Е * ЗГ2"

_

и1 + 121а дС4

4 д4и1 2 д2и1 ^ , ди1

а

дС 2 +)+ ^ —

— V! + /

Е *Н*

2vw0 + (1 — V )а

Р/ —

ди0

5С"

2 ди0 ди1 а2 / ди0

— + Т

Д0 0.5 — 1п

1Рдт ей*

4^

д2и1

"Е* дт2"'

Текущие координаты осевой линии выражаются через функции (8):

, ч , ч „ еД0 ^0 ч

ж( — ,£) = Ж0(з) + — Д0«0 + (V1 + и1),

/ ,ч / \ , „ еД0 _ ч

у(М) = У0(—) + "^—^0^0--2—1 + ад1),

*■'> = 5 (0

-1

и1

^ дv1 + ди1

1 2 Vдт 1 дт Однородные начальные и краевые условия преобразуются к виду

ди0 ди1

и0 = и0 = v1 = и1 = 0,

дС

ди1

и1 = v1 = и1 = 0, -7— = -7—

дv1

37

дС "

ди1

0 при £ = 0, С = ^ и любом т,

(9)

(10)

0 при т = 0 и любом (11)

дт дт дт 1

Уравнения (6), (9) с начально-краевыми условиями (11) и соотношениями (10) представляют собой приведенную к одномерному виду математическую модель движения трубопровода как оболочки при конечных деформациях.

2

1

В качестве конечного результата нас интересует перемещение осевой линии трубопровода вдоль нормали к ее начальному положению: = п0), где w — вектор перемещения, а п0 — нормаль, которые в системе отсчета хОу выразятся так:

, ч. , , ч. йуо. ^Хо.

w = (х - Хо)1 + (у - УоЛ, по = 1 -

смысл х, у, х0, у0 пояснен выше. Отсюда, учитывая свойство естественной параметризации кривой х02 + у02 = 1, запишем

М = ^ (^1 + ^1).

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

м-п = — ^ (VI + ^1). (12)

4. Деформации стенки трубы и критерий прочности

При выводе уравнений равновесия оболочки в [11] деформации трехмерного упругого тела представлены в виде

« -1 + Ка(К — Ко), -вв ~ -2 + К2(К — Ко), £sв « и + т(К — Ко).

Здесь величины К — Ко имеют порядок толщины стенки Н. Нас интересуют главные порядки [11]:

1 du 1 dA 1

e1 = т^;—+ ~п=; + ki w —

Ads АБдв 1 2A2

/ dw\ / dv V V Ids) + V ds)

1 ОВ 1 ОУ ,

-2 = АВа7и + В~т + к2Ш,

Ад (и \ ВО (V \ ..

и = Вт (а) + А¥Лв)- (13)

Учитывая значения параметров поверхности оболочки (2), переходя к безразмерным переменным из (3) и раскладывая в ряд по -, получим

du0 а2 /dwo\2 . Л Г du1 2dw0 dw1 ( 2(dw0\2 с

а—— а -тгг^гг- + / w0 + а —— — а-

v.^0 а dwo , , . „

e1 = — TV l^1 + e sin 0

dZ dZ dZ 4 \dU dZ

e2 = w0 + e sin6(w1 — v1), w = e cos 0 ^u1 + адv1" — /u0^ • (14)

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

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

е* + £у ±\ (е* - £у)2 + тХу < 2[е],

(15)

обозначения обычные. По аналогии для трубы можно ввести критерий прочности, заменяя в (15) ех ^ в\, еу ^ е2, 7*у ^ ш. Тогда получим

е1 + е2 ±у/(е1 - е2)2 + ш2 < 2[е].

(16)

Это соотношение имеет смысл условия максимальных деформаций удлинения. Совместно с формулами (14) условие (16) может при принятых предположениях относительно поведения механической системы применяться для оценки прочности стенки трубопровода.

5. Алгоритм нахождения численного решения

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

Для решения системы уравнений нулевого приближения (6) использован пакет МаЛеша^са 4.2, разработка программы решения методом последовательных приближений не представляет больших затруднений [15].

В качестве численного метода решения системы уравнений первого приближения (9) выбрана явная трехслойная разностная схема с уточнением значений, способ построения которой подробно изложен в [16]. Суть уточнений в том, что вычисляется по строго явной схеме, при определении ь используются уже вычисленные значения на верхнем слое, и1 вычисляется по неявной схеме методом прогонки.

Разностная схема по пространственным переменным строится известным интегроин-терполяционным методом [17]. В результате этого построения для средних точек интервала изменения переменной £ получаются следующие расчетные формулы:

7+1 = 2^' - V'-1 +

Е*а2 й2 1-

V

РЛ2ш2 й2 2

«1 -+ -

Е *й22

Р^ш2

7

V —

2^22

и,

*з

й*Я3Р*ш2

0,5 - 1п

!Рдт еи/ Д

4^

Е * й2 Г 1 + V а , , , Л ,

+ рдшл—^ - и*-1) + ь+

2

Ь0г - ; (и0(^+1) - и0(г-1}) - 4й2 (Ь0(г+1) - ь0(г—1)) К+1 - ^1)

"С

ь:

1,7+1 = 2ад7 — 1 —

Е*й2 ( 7 й*2 Га4

р4Я2ш2 1

Ь + 12

( Ь+2 - 4—+1 + - 4ь^-1 + Ь-2 I -

г—2

а

- 1 ьг+1- 2ьг + 1

+ (и7+1- и7—^ - + ^

С

2vWог +

,

а

+ (1 - ^

а2 (

«0(1+1) — М0(г-1)

ЩУ '0('+1) — '0('-1)

и

Е *к*

а

4к2 ^ '0(т)

Р/г У/ —

'0('-1) '+1 — '¿-1 +

2^и

Д0 0, 5 - 1п

^Рдт еи{ Д

4^

2к

гз + - г, + '3+ - ад, .

После этих расчетов можно вычислить и\+1 по неявной разностной схеме

(17)

,+1

и

Е *а2

Н2

РгЩш2 к2

иг-+1 2и3 + иг-+1

+

1 — V Е *Н2Т ~ РьЩ^2

и

3+1 = 2и, — и,-1+

Е *к2т ргК0)и2

1 + и а ( ,+1 ¿+1\ , иа , _ с -+1 — ^ + 2кс ^

4 к

— Л') +

1 — V а2 ( \ а ,

и0' — 2кг (^0('+1) — 2и0г + ^(¿-1)^ + (1 — V) — ('0('+1) — '0('-1)

а

2к3

ад0('+1) — 2woi + ад0(г-1^ (,11 — ад'-Л + ( ^(¿+1) — '0('-1)) ( ад'+1 — + ад'-/

,+1

,+1

,+1

,+1

+

За3 2к|

/ ('0('+1) — '0('_1)) ('0('+1) — 2'0' +

(18)

Здесь кт — шаг по времени; к^ — шаг по координате; = гк^; /з = /((¿, т,); и/ = т,);

Щ' = и0(('); '0' = '0«'); из = Й1((',т,); гз = ^(С',^); ' = '1(6,^). Расчетные формулы на краях интервала при г = 1, N — 1 немного отличаются от (17), но отличие незначительно, а метод учета краевых условий хорошо изложен в [17], поэтому здесь они не приводятся. Система (18) решалась методом прогонки, учет однородных краевых условий

3+1

также прост: и0

= 0.

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

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

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

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

— блок ввода исходных данных и задания разбиения по координате £;

2

1

з

— блок расчета текущего шага по времени;

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

Окончательные результаты формировались в виде текстовых файлов, которые затем обрабатывались пакетом МаЛеша^са 4.2 для графического отображения.

6. Результаты тестовых расчетов

В качестве физических и геометрических параметров модели выбраны следующие значения величин: ра = 105 Па, рдг = 1700кг/м3, V = 0.31, Е = 6.8 ■ 1010 Па, р4 = 2700кг/м3, к0 = 5 м, к = 0.005 м, Е0 = 0.3 м, Ь = 6000 м, р7 = 850кг/м3, в = 27.38кг/м4, ^ = 2000 Па-с, уз0 = 2 м/с. Эти параметры соответствуют потоку нефти в алюминиевой трубе, погруженной в мягкую почву.

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

х меняется от 0 до 6000 м. Этот пример является тестовым и выбран для проверки математической модели.

В качестве конечных результатов выводились численные значения перемещений осевой линии (12) как функции £, т и значения сдвиговой деформации ш (14) как функции £, 0 в конечный момент времени. В выбранном примере значения деформаций продольного и окружного растяжения е1, е2 оказались почти не зависящими от координаты 0.

График зависимости и-п = приведен на рис. 3, рассчитан интервал времени

семь часов. Малый интервал времени расчета вызван сложностью модели и использованием персонального компьютера РеП;шш-4 1700 МГц. Устойчивый счет достигнут при шаге по времени 0.000085 с.

График зависимости ш = ш(з,0) приведен на рис. 4. На нем видно, что наибольшие искажения углов при выбранной форме осевой линии возникают в окрестности того конца

у = 10-8 ■ х(х - 3000)(х - 6000),

У 150

и

100

50

0

х

-50

-100

-150

0

1000 2000 3000 4000 5000 6000

Рис. 2. Форма осевой линии в начальной конфигурации.

Рис. 3. Поперечное перемещение осевой ли- Рис. 4. Деформация сдвига как функция в, в. нии.

Рис. 5. Изолинии деформации сдвига.

трубопровода, в сторону которого движется жидкость. На рис. 5 показан тот же график в виде изолиний и = const.

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

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

[1] Жуковский Н.Е. О гидравлическом ударе в водопроводных трубах. М.; Л.: Гос. изд.-во техн.-теорет. лит., 1949. 104 с.

[2] Члрный И.А. Неустановившееся движение реальной жидкости в трубах. М.: Недра, 1975. 296 с.

[3] Вольмир А.С. Оболочки в потоке жидкости и газа: задачи гидроупругости. М.: Наука, 1979. 320 с.

[4] Berger S.A., Talbot L., Yao L.S. Flow in curved pipes //Ann. Rev. Fluid Mech. 1983. Vol. 15. P. 461-512.

[5] Клочков Б.Н., Кузнецова Е.А. Нелинейные режимы изменения формы упругой трубки с потоком жидкости в ней // Механика жидкости и газа. 2000. № 4. C. 46-55.

[6] Lynch D.G., Waters S.L., Pedley T.J. Flow in a tube with non-uniform, time-dependent curvature: governing equations and simple examples //J. Fluid Mech. 1996. Vol. 323. P. 237-265.

[7] Lee U., Pak C.H., Hong S.C. The dynamics of a piping system with internal unsteady flow // J. Sound and Vibr. 1995. Vol. 180, N 2. P. 297-311.

[8] Ткаченко О.П. Кинематика и динамика подземного трубопровода при конечных перемещениях // Вычисл. технологии. 2003. Т. 8, № 4. С. 97-107.

[9] Рукавишников В.А., Ткаченко О.П. Численное и асимптотическое решение уравнений распространения гидроупругих колебаний в изогнутом трубопроводе // Прикладная механика и техническая физика. 2000. Т. 41, № 6. С. 161-169.

[10] Лурье А.И. Нелинейная теория упругости. М.: Наука, 1980. 512 с.

[11] Власов В.З. Общая теория оболочек и ее приложения в технике. Избр. труды. Т. 1. М.: Изд-во АН СССР, 1962 (С. 15-439).

[12] Лойцянский Л.Г. Механика жидкости и газа. М.: Наука, 1987. 840 с.

[13] Феодосьев В.И. Сопротивление материалов. М.: Наука, 1986. 512 с.

[14] РжАницын А.Р. Строительная механика. М.: Высш. шк., 1991. 439 с.

[15] Дьяконов В.П. Математика 4: Учеб. курс. СПб.: Питер, 2001. 656 с.

[16] Ткаченко О.П. Построение математической модели распространения волны давления в изогнутом подземном трубопроводе и приближенное решение ее уравнений: Дис. ... канд. физ.-мат. наук. Владивосток, 1998. 141 с.

[17] Самарский А.А., Гулин А.В. Численные методы. М.: Наука, 1989. 432 с.

[18] Власов В.З. Тонкостенные упругие стержни. М.: Физматлит, 1959. 568 c.

Поступила в редакцию 15 июля 2005 г., в переработанном виде — 25 октября 2005 г.