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

ч> На правах рукописи

г

Казанцева Евгения Владимировна

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

05.13.16 - применение вычислительной техники, математического моделирования и математических методов в научных исследованиях

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук

Барнаул - 1997

Диссертация выполнена в Сибирской государственной геодезической академии (Новосибирск)

Научный руководитель: доктор физико-математических наук,

профессор В.Н.Врагов

Официальные оппоненты: доктор физико-математических наук,

профессор А.Ф.Воеводин

Ведущая организация: Вычислительный центр СО РАН, г.Новосибирск.

на заседании диссертационного совета Д 064.45.02 в Алтайском государственном университете (6560ЭЭ, Алтайский край, г.Барнаул, ул. Димитрова 66).

С диссертацией можно ознакомиться в библиотеке Алтайского государственного университета.

кандидат физико-математических наук, О.Б.Бочаров

Защита диссертации состоится "ЛУ 1997г. в {£_ часов

Автореферат разослан

Ученый секретарь диссертационного совета доктор физ.-мат. наук, профессор

С.А. Безносюк

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

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

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

На защиту выносятся следующие результаты:

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

2. Численный алгоритм для реализации предложенной модели.

3. Аналитические решения: задачи об образовании полусферической и полуэллипсоидальной газовой полости; задачи об обтекании изменяющегося во времени эллипсоида вращения.

4. Результаты численных расчетов, полученных по предложенной модели и численному алгоритму.

Научная новизна:

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

- 3 -

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

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

4. Получено аналитическое решение задачи при полусферическом и полуэллипсоидальном изменении газовой полости. Выполнен анализ полученных решений.

5. На основании решения конкретных задач и сопоставления результатов численного решения с экспериментальными данными, известными из литературы, показаны: характерные особенности возникающего течения жидкости и движения газа в полости; надежность предлагаемой теоретической схемы и численного алгоритма.

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

Апробация работы. Основные положения и результаты работы докладывались и обсуждались на XXVI Всесоюзной научной студенческой конференции (Новосибирск, 1988); на II Конференции молодых ученых Сибири и Дальнего Востока (Новосибирск,1989); на конференции по проблемам создания бессточных систем водоснабжения промышленных предприятий и экологии Средней Азии (Ташкент, 1991); на научных семинарах кафедры теоретической механики Новосибирского гос. университета; на научных семинарах по проблемам фильтрации в ИГиЛ СО РАН; на научных семинарах по проблемам механики сплошных сред в ИГиЛ СО РАН; на научном семинаре МВТ СО РАН (Новосибирск, 1994);на 51-ой научно-технической конференции Новосибирской государственной академии строительства (Новосибирск,

_ 4 _

1994); на научно-методическом семинаре Сибирской государственной геодезической академии (Новосибирск,1994); на Международном симпозиуме по гидродинамике судна (С.-Петербург, 1995); на научном семинаре ВЦ СО РАН (Новосибирск, 1996).

Публикации. По теме диссертации опубликовано 5 работ. Объем работы. Диссертация состоит из введения, четырех глав и выводов. Содержит 148 страницы, включает 39 рисунков и 6 таблиц. Список литературы включает 71 наименование.

СОДЕРЖАНИЕ РАБОТЫ Во введении дается литературный обзор по теме исследований, обосновывается новизна и актуальность рассматриваемой проблемы. Еормулируется задача исследования и излагается краткое содержание работы.

В первой главе диссертации дано описание физико-механических процессов, возникающих при истечении газа из донного отверстия (дно-горизонтальная плоскость) в несжимаемую жидкость с образованием в начальной стадии процесса газовой полости. Методом осреднения по площади "живого" сечения (рисЛ) трехмерных уравнений газовой динамики получена система квазилинейных дифференциальных уравнений, описывающая неустановившееся течение газа в вытянутой по вертикали полости, имеющей переменную шю-дадь сечения А=1Ш2(г^). Движение жидкости, вызванное возникновением полости, предполагается безвихревым и описывается уравнением Лапласа, которое записывается в цилиндрической системе координат. Решение этого уравнения ищется в области П2, причем эдна из граничных поверхностей (поверхность раздела между газом ,1 жидкостью) меняется и по времени, и по пространственным ■соординатам.

Математически задача сформулирована так: найти решение уравнения Лапласа

в области 02(см.рис.2) при граничных условиях:

- 7+ наТ4:ф=о,

на т2: на §§=0, (2)

на д| + '

Ш+ + ^ - ё(н-2)=0 14311 г-н(в-*)'

ря

(3)

при г=Н(гД); и системы уравнений, описывающей движение газа: ^(р.,А) + д^р^А) = (4)

и2

^(Р1А(е + -4-))+^(р1^1А(в + + (6)

6=0^, р^р^Т.,, (7)

в области П., :{СК2;^к(-Ь);0<Шк} с гранич1шми условиями: на левой границе области при 2=0

0=^(1;), т1=*3(±) (если |и,|^с); (8)

о=14(1;), (если |^|<с); (9)

на правой (подвижной) границе области при й=ак(1;) Да,

с начальными условиями для искомых функций: ф(г,2,0)=фо(г,г;), ,

, , (II)

Н(г,0)=Н°(г) , Т1 (2,0)=а?°(г) .

Таким образом, математическая задача об определении параметров газа и жидкости при выбросе газа в несжимаемую жидкость свелась к нахождению функций ф=ц>(г,%,Ь), г=Я(гЛ), г=Н(г,1;), р^р^аД), и^тл, (я,1;), (г,1;),удовлетворящ0Х уравнениям (1),(4-7)» граничным условиям (2),(3),(8-10) и начальным условиям (II).

Здесь и далее по тексту использованы следующие обозначения: 1;-время; г-ось цилиндрической системы координат, направленная вертикально вверх; г-радиальная координата; ср=ф(г,г,1;) - потенциал вектора скорости жидкости и2; г^- средняя по сечению скорость газа; р^-плотность (1=1 - газ; 1=2 - жидкость); р^-давление; ра-атмосферное давление;е=с^т.,- внутренняя энергия газа; су- удельная теплоемкость газа; и - универсальная газовая постоянная;

о

¿^(гД)- площадь поперечного сечения полости; г=БЦг,1;) -поверхность раздела между газом и жидкостью; г=н(г^) -ордината свободной поверхности кидкости; с=р1 г^ А - массовый расход газа; с=У7й -скорость звука в газе; 7-показатель

о '

адиабата; ^(-Ь)- траектория вершины полости; £¿(1;)- заданные Функции; ^ускорение силы тяжести.

Граничные условия для . потенциала <р есть условие непротекания на твердой стенке, условие симметричности потока, кинематические условия на свободных поверхностях. Условия (3) выписаны с учетом равенства давлений на поверхности раздела двух сред и интеграла Кош-Лагранжа.

При выводе уравнений (4)-(6) для определения средних по сечению величин применена формула

~ ДКх.у.г.ЮаА.

А

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

^ ЗА _ „ oiroÖR

учитывающего работу силы давления на расширение (сжатие) газовой полости.

Установлена принадлежность системы дифференциальных уравнений (4)-(6) к гиперболическому типу с характеристиками

число условий (8)-(10)' для уравнений (4)-(6) на границах (z=0; z=zk(t)) определено, как это принято для гиперболических уравнений, числом "входящих" в область характеристик.

Здесь, по-видимому, следует отметить, что А=0 при z=zk(t) (вершина полости) не вызывает особых затруднений, т.к. возникающая особенность в уравнении (4) (или (6)) д п + 9 ii ^ - 2p1 fÖR + ■>. ÖR1

9tPi + 9z(Pivi > ~ ~ "TrlcFf + vim является устранимой. Действительно, с учетом условий (2), (10) и фг= | г>2| sine (6-угол между вектором üj2 и осью 0Z), при R-0 ¿в -О, получим

1 fÖR ÖKl^r" .

К + Щ= * 3

R-»0 &Ч) . K ' R-0 pR p zk

Как видно из граничных условий (8)-(IQ), в газовой

- а -

полости теоретически возможны любые (сверхзвуковые, дозвуковые и переходные) режимы течения. Отсюда и возникает вопрос о качественном поведении интегральных кривых для полостей, форта поверхности которых известна по данным физических экспериментов (С.С.Кутателадзе, В.Е.Накорякоз, КЕМ г.Миасс). Ответить на поставленный вопрос в полном объеме - значит решить в общем виде уравнения газовой динамики. Наш был рассмотрен случай движения полости с постоянной конечной скоростью и. Переходом к переменной. £=г-тп, получены обыкновенные дифференциальные уравнения, для которых построено поле интегральных кривых. На рис.3 представлены распределения числа Маха М=М(£) для г=к(£,-Ь). Как видно из данного рисунка, существуют решения, соответствующие сверхзвуковому (кривые (I)), дозвуковому (кривые (2)) течениям. Две кривые (3),(4), проходящие через точку сужения полости, обеспечивают плавный переход от сверхзвукового (дозвукового) к дозвуковому (сверхзвуковому) режимам течения. Имеются решения, являющиеся чисто математическими (замкнутые и возвратные кривые (5)).

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

и работах С.К.Годунова и др., Б.Л.Рождественского, Н.Н.Янен-ко, Г.И.Марчука, А.А.Самарского, Р.Р.Рихтмайера, Дк.Мортона, А.Ф.Воеводина, И.К.Яушева. -

Для аппроксимации системы уравнений (4)-(6) строится "подвижная" разностная сетка с шагами 1хп+1 =2^.((п+1 )т)ЛЛ и т по координатам г и t соответственно. Во внутренних узлах сетки уравнения аппроксимируются четырехточечной неявной итерационной разностной схемой ф^1- ф« г1*1- ф^,1- ф™,1

-г т А»+1 А**1 з ' 1

0=1,.., 1^-1 ,

Р., А Р-, V* 0

ф = р1А(е+^) ; <|> = р^А + р., А А Атл, (р1 (е+21-)+р1) ; ? = гч <ЭА р1 ш „ дА

К разностным уравнениям (13) присоединяются разностные аналоги граничных условий (8),(9),(Ю) и разностные уравнения, аппроксимирующие схемой "уголок" уравнения, записанные в характеристическом виде. При этом берутся те уравнения, которые соответствуют характеристическим направлениям, "приходящим" к границам изнутри области.

Разностные уравнения (13) решаются методом матричной прогонки после их предварительной' линеаризации по формуле

фГн-1 фП+1 .V + стп+' ___

Щ г ] •

Й={р1,У1,Т1}.

При численном интегрировании уравнения Лапласа,

описывающего движение жидкости в области с целью

^спрямления" границы, вводятся новые независимые переменные

Для аппроксимации уравнения Лапласа,

записанного в новых переменных

32ф Э2ф б2ф бф

—^ - 2Ъ- + С—« - <1— = О, (14)

ЗГ д^дг] дтГ 0т\

ге 1+гр 1 1

где с=—2^-, - 2Ьг?Т) + сг^ - р, используется

двухшаговая схема стабилизирующей поправки на неравномерной сетке. Сетка сгущается в точках замыкания границы полости с осями £,т}. Разностная схема является схемой, полной аппроксимации. На каждом шаге разностные уравнения решаются методом прогонки.

Вычислительный процесс, при переходе с п-го временного слоя на (п+1 )-й, организуется так: вычисляется положение границы раздела, определяются параметры газа, методом установления решается уравнение Лапласа. На каждом временном шаге указанный процесс повторяется (итерируется) до достижения заданной точности вычислений:

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

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

Изменение газовой полости, • имеющей форму полусферы,

рассмотрено в первом разделе третьей главы. В случае

сферической симметрии для потенциала ср имеем: 2 •

, где а2=Ц2+г2, р2=г2+г2 (15)

а=а(-Ь) -радиус сферы, ¿=¿(1;)= ^ -скорость движения границы раздела.

Т зт

Давление газа в полости: Р-1=Ро(йг) .здесь

7-показатель адиабаты газа, т=т(1;)-масса вдуваемого газа, • т0-масса газа в пузыре в начальный момент времени, о0-радиус пузыря в начальный момент времени. Из интеграла Коши-Лагранжа получено уравнение, описывающее расширение сферической полости при заданном изменении массы газа:

- о.

Откуда при 7=4/3 и при га=т0 (расширение полости определяется начальным вдувом газа) получена формула, дающая связь между временем и радиусом сферы:

р- = (1 - |п - 1п2)/гп"\

О -3 5

ЗДеСЬ С2=р0/р2, П=а/ао-1.

Ео втором разделе решена задача для случая, когда полость изменяется по времени, сохраняя при этом форму полуэллипсоида вращения. Для решения задачи сделан переход к эллипсоидальной системе координат (Х,ц), которая связана с цил-линдрической (г,г) по формулам: г=сЛц, г=су(Х2-1 )(1-ц2), где

с2=а2-ь2,а а=сг(1;) И Ь=Ь(-|;) -ПОЛУОСИ ЭЛЛИПСОИДЭ. СфОрмуЛИрОвана постановка задачи в переменных X, р.. Для потенциала ср получено выражение:

"п-^+С ба -4С Ь

с о—с

3ь

Для проверки правильности полученного решения рассмотрен предельный случай: при ь-*а, на основании следующих предельных соотношений:

с-О, сА.-Н=Уг2+г2, 1п^-=2^ + + ...], при ь-к7, показано, что из формулы (17) вытекает формула (15). Вычислены распределения давления в жидкости и по поверхности тела Р2

Р2

и ,, ^о.УУ ДЕ, 2 1 ^ 1 2_Й2 >..2.1/2]

Здесь 02=1(ЗА.2-1 - - ' полинош

Лежандра 2-го рода; А,0=^ - поверхность эллипсоида в координатах А,, ц.;

СЮд

т? -ЯН. V -551 V -сь£ • В - О

о о 2 .. ..0 .. -0

Разность давлений в вершинах эллипсоида :

х=х

£2 р2

р.

,р2 0,(А.0) «ЛЬ,

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

Для функции тока и кинетической энергии жидкости получены следующие выражения:

♦ ■ - ье

В = «ргсьЗ(0о(ло,'4,§й£; ♦ Ц!. ъ

На основании формул предельного перехода, показано, что

Е-*27Ср2а3а2, (Ь-а).

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

С1п^с--^Г 2 а_с ь2-»

Построены линии (р=сог^ и <{>=оопб1 (гидродинамическая сетка) при скорости набегающего потока иоо=4м/с и ог/ао=1+21;2, ь/ьо=1+|;2, «о=-ьо. На рис.4 представлены соответствующие моментам времени t=0.5c; t=I.5c; t=2.5c расчетные распределения ф и ф (рисунки в верхнем ряду), и линии р=оош^ (нижний ряд). Распределение давления на поверхности эллипсоида и на оси потока представлено на рис.5. Здесь сплошная жирная линия соответствует случаю набегания потока на расширяющийся эллипсоид; тонкая сплошная линия -случаю расширения эллипсоида без набегания потока (0^=0); пунктирная линия -случаю обтекания эллипсоида без его изменения. Эти графики наглядно показывают неравномерность распределения

давления на поверхности эллипсоида.

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

Для проверки точности предложенных численных алгоритмов все три программных модуля были протестированы автономно и в комплексе на известных аналитических решениях. Для уравнений газовой динамики в качестве тестов были взяты: задача о распаде произвольного разрыва (результаты представлены на рис»6); задача о течении газа в сопле Лаваля и др.

Для проверки точности решения уравнения Лапласа были использованы полученные в данной работе аналитические решения о полусферическом и полуэллипсоидальном изменении газовой полости. Результаты решения задачи о полуэллипсоидальном изменении полости представлены на рис.7, на которых изображены поле давления (рисунки в нижнем ряду) и поле течения жидкости - линии ср=сог^ и ф=сог^ (верхний ряд) для различных моментов времени. Полученное по численному алгоритму значение потенциала ф отличается от вычисленных по формуле аналитического решения не более- чем в пятом знаке после запятой. Далее, была выполнена серия расчетов для полной задачи (1)-(П). Результаты, полученные при этом, сравнивались с экспериментальными данными. Геометрические размеры

пузыря (см.рис.8), полученные численным методом, отличаются от полученных методом фотовизуализации приблизительно *На 7-9%, что свидетельствует о достаточной точности математической модели и численного алгоритма.

Основные результаты диссертационной работы:

1. Предложена математическая модель для описания процесса подводного выброса газа из донного отверстия в несжимаемую жидкость с образованием вытянутой по вертикали газовой полости. Данная модель включает уравнение ■ Лапласа, уравнения неустановившегося течения газа в полости, необходимые граничные и начальные условия. Модель предназначена для определения: параметров неустановившегося течения жидкости, вызванного возникновением полости, поверхности раздела, распределения параметров газа по времени и по длине полости.

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

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

4. Получены формулы для определения поля течения и поля

давления при обтекании пульсирующего эллипсоида потоком несжимаемой жидкости.

5. Выполнены расчеты, показывающие: точность численного алгоритма характерные особенности течений жидкости и газа при подводном выбросе газа; применимость математической модели для решения конкретных технических задач.

Основное содержание работы опубликовано в следующих • работах:

I.Чернышева Е.В. Динамика кавитационной газовой полости с учетом волновых процессов в газе//Материалы XXVI Всесоюзной студенческой конференции "Студент и научно-технический прогресс", Новосибирск, апрель 1988г.:-Новосибирск, Изд-во НГУ, 1988. -С.73-79.

2.Чернышева Е.В. Продольное обтекание изменяющегося во времени эллипсоида вращения в безграничной жидкости// Труды II конференции молодых ученых Сибири и Дальнего Востока, Новосибирск ,1988г.:-Новосибирск,Изд-во НГУ,I988. -С.180-184.

3.Тимеркаев P.M..Чернышева Е.В. Моделирование фильтрации подземной жидкости в зоне влияния накопителей промышленных стоков//Межвузовский сборник научных трудов "Математический анализ и дискретная математика". - Новосибирск: -Йзд-ео НГУ, 1989. -С.154-170.

4.Тимеркаев P.M..Чернышева Е.В..Чернышева Р.Т. Задачи организации защиты подземных вод от загрязняющего воздействия накопителей промышленных стоков// Труды конференции по проблемам создания бессточных систем водоснабжения промышленных предприятий и экологии в регионе Средней Азии, Ташкент, май 1991г. - Ташкент, 1991. -С.82-84.

5-Врагов В.Н. .Чернышева Е.В. Продольное обтекание перемещающегося во времени эллипсоида вращения потоком потенциальной несжимаемой жидкости//Труды Международного симпозиума по гидродинамике судна, С.-Петербург, май 1995г.:-С.-Петербург, Изд-во НЛП "Форма", 1995. -С.289-304.

Поперечный разроз газовой полости

у

2Л>

Л,

Х=-К„(2,£) \ ) х=Ыг,{.) '

Рис. 4.

Обметь решения задами

V Х-

V N

V N Ч N V ч Ч ч, N г, 2 = Н(ХД)

- Л, N г,

- ------ .. . -РИС. 2. • %

Качественное поведение интегральных кривых в газовой полости

• *

Ч

?

рис.3

«и -» о а ч • « ю

Линии с^сопог и ф=сопз* (верхний ряд) и линии.равных давлений-при обтекании расширяющегося эллипсоида.

рис.4

—— при набегании потока на расширяющийся эллипсоид; - при расширении эллипсоида (1^=0);

---- при набегании потока на покоящийся эллипсоид.

рис.5

П Г-