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

ББК В 253.322 УДК 532.546

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ РАСПРОСТРАНЕНИЯ ИМПУЛЬСА ДАВЛЕНИЯ В УПЛОТНЯЮЩЕЙСЯ ПОРИСТОЙ СРЕДЕ.

Лукин С.В., Урманчеев С.Ф.

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

Введение

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

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

1. Математическая модель

Рассмотрим дисперсную систему, состоящую из плотно упакованной зернистой среды, поровое пространство которой заполнено жидкостью. Для описания распространения волн в таких средах примем допущения аналогичные [2], [3].

В рамках принятых предположений для изучения распространения волн в таких средах используем двухскоростную модель насыщенной пористой среды [2], [3]. Система уравнений модели пористой среды, насыщенной жидкостью вместе с замыкающими соотношениями имеет вид [2], [3]:

Ф1+Ф^ = 0 Ф 2 , Ф 2^2 = 0

д? дх ; д? дх ;

др, ^ й9у9 дрЛ ^ да 2-

р1“^ = ^1"^“ ^2 Р 2^г _ 2 ^12 +-^-

а? дх ; а? дх дх (1)

й2а 2* 17 Ш 2* _ 1 („ 17 о ) й2£ дУ2

—Т---------КГ~Г~ _ “~(а 2* “ Ье*° 2* ) ^

а? а? ?20 . ш дх ;

а1 2 _ 1 ; р1 _а1р1 ; Р 2 _а 2р 2 ;

^ Р2*

р1 _ р0 + ^1 (р1 “ р10 ); р2 _ р0 + С2 ( 2 “ Р 20 ); а 2 , где р2* _ “Уа 2* .

р 0

Здесь р I , р I, а 1, р,, - истинная и приведенная плотности смеси, объёмное содержание, давление и

о а

скорости фаз соответственно; &2* - продольный компонент деформации твёрдой фазы; 2* - приведенное

напряжение в скелете; 1 _ 1 относятся к жидкости, 1 _ 2 - к твердой фазе. Модули упругости скелета определим через скорость звука в скелете, соответствующими мгновенному и длительному модулям:

*

Лукин Сергей Владимирович - аспирант Института Механики УНЦ РАН

Урманчеев Саид Фёдорович - канд. физ. - мат. наук, доцент кафедры механики сплошных сред БашГУ

И,

где 1 , «

Ef * _ р 20* Ее, _ р 20Ое,

Ег * Е„,

(2)

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

статический модули упругости скелета пористой среды; индекс нуль соответствует начальному состоянию.

р

Силу межфазного взаимодействия Р12 представим в виде суммы силы вязкого трения Стокса н и силы

т7

присоединённых масс ™ [2]:

Е _ Е + Е

1 12 н г”

где

а1ч1 а 2ч2

а? а?

(3)

где Н1 - динамическая вязкость жидкости; а2 - радиус частиц, составляющих скелет пористой среды.

Коэффициенты Лн и характеризуют вязкое и инерционное взаимодействие фаз и зависят от структуры

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

режимы с ^е12 ~ 10 и коэффициент Л н можно полагать постоянным.

Исследуем эволюцию волнового процесса с монотонно убывающей пористостью по закону:

(4)

а1 _ К • х +М • х+М

где К = -0,29,М = 0,291, N = 0,627. Здесь пористость убывает от 0,7 до 0,3.

2. Конечно-разностная аппроксимация системы урав нений насыщенных пористых сред

При заданных значениях л ™*, Ие*, ?20 выписанная система уравнений (1) является замкнутой. Процедура, используемая для её численного интегрирования, на каждом шаге по времени состоит их трёх этапов. На первом этапе «замораживается» межфазное силовое и тепловое взаимодействие, т.е. полагается М = 0, и дифференциальные уравнения системы интегрируются с помощью двухшаговой схемы Лакса-Вендроффа

[5].

На первом шаге вычисляются значения переменных для промежуточного (п+1/2)-го временного слоя:

т , т

п+1 1 X

и 2 _1 (и ”

.1 о ^ У

2 ' ,где] = 1,2,$,N-1-, (5)

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

]+1

и3 -

) А?

2Дх

- ¥М ¥М ¥М fft fft ¥М —.

^+1 - К, + 2 (г/+1 + Vя (,:■+1 - к)

3

3

А?

Ах

1

1

V 12

V ■/+2

з-

1 1 1

пл— пл—

w 21 - w 2

,, з+— ,, з-

2

2

V

2

где] =2,3, $ ,N-1. (6)

На втором этапе производится учёт членов межфазного взаимодействия по схеме первого порядка точности:

и "+1 _ ип

-Д?0\и

(7)

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

фГ _(1 - * >Г + ч № +*;+), (В)

(ч )"+1 - 2 (ч Г +(’■ К\> Н (ч )”1 - (ч )Г|&| (ч)" - (ч)

-в противном случае (9)

Ч _

п+\

I.

> 0

У = 3,4, ... , N-2,

Ф = {р 1, Р 2’ Рі^1’ Р 2^2’° 2*’£ 2 }, где У = 3,4, ... , N-1. (10)

Шаг по времени в численном алгоритме выбирается из условия устойчивости решения и точности расчётов процессов межфазного взаимодействия:

А? < тт(д?™, А?[ )

і=1,2,

(11)

А?крь = -

Ах

С,

а?: = ■

Л ^1

- время силовой (стоксовой)

где I '11 ' - условие Куранта-Фридрихса-Леви;

релаксации при межфазном взаимодействии; Ах - Шаг по координате.

Для задания граничных условий введём фиктивные ячейки, соответствующие узлам у=1 и у =М. Границы физической расчётной области проходят: левая - между первым и вторым узлами, а правая - между N-1 и N узлами. Использовалось условие протекания на левой и правой границах:

ФГ = Ф 2И Ф = {Р1’ Р 2’ РЛ’ Р 2У2’° 2*’£ 2 }.

(12)

На контактной границе жидкость - пористая среда (О) задаётся скачок параметров, характеризующих

0

^ О V Р

скелет, - объемного содержания твердой фазы 2, ее плотности У 2 , скорости 2, деформации 2* и

приведенного напряжения ® 2* . Предполагается, что перечисленные параметры слева от границы равны нулю. Кроме того, должно выполняться условие равенства полного напряжения слева и справа от этой границы

Рю- Рю+ °

2*в +•

(13)

Очевидно, что должно быть справедливо равенство ® 2*0+ _ , которое представляет собой условие

Терцаги и откуда следует, что (13) выражает условие непрерывности давления в жидкости.

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

3. Численные результаты

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

Перейдём к рассмотрению численных решений, полученных с помощью вышеописанной методики. Численные исследования распространения импульса давления в пористой среде насыщенной жидкостью проводились в соответствии со следующей схемой. Ударная труба длиной около полутора метров, примерно на две трети заполненная кварцевым песком, доверху заливалась водой. Свободная поверхность жидкости нагружалась импульсом давления. Средний радиус песчинок в расчетах принимался равным а2 = 0,25 мм,

а1П = 0.3

* = 1700 м/с

а значения модельных параметров і ;

В* = 850 м/с

пористость среды

^20 _ ^и0 Мкс , Лц 6°°, _ °.75 . 2 соответствии с указанной экспериментальной схемой была построена

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

Рассмотрим волновое течение в насыщенной пористой среде при переменной пористости, когда она непрерывно уменьшается от 0.7 до 0.3.

Рис. 2. Изменение максимального значения безразмерного давления в жидкости в «быстрой» (верхний график) и «медленной» (нижний график) волн при убывающей пористости (сплошная линия) и при постоянной пористости

(пунктирная линия, /20 =690 мкс). Значения параметров: ® 1 =0,7 убывает до 0.3, Л ^ = 600, Л т =0.75.

На Рис. 1 представлены эпюры давления при прохождении импульса из жидкости в насыщенную

жидкостью пористую среду с пористостью, убывающей от значения ^1 = 0,7 на границе раздела до ^1 = 0,3 в конце расчетной области. Числовые указатели у кривых соответствуют времени в мс.

В отличае от аналогичного случая при постоянной пористости (а1 _ 0,7), при уменьшении пористости, амплитуды давления жидкости в «быстрой» и в «медленной» составляющих ударной волны заметно увеличиваются (за время от 0,6 мс до 0,9 мс амплитуда быстрой волны увеличивается в полтора раза).

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

Случай изменения амплитуды давления жидкости в «быстрой» и «медленной» составляющих волны представлен на Рис. 2. В «быстрой» волне при упругом скелете и вязкоупругом с 1;20 = 690 мкс амплитуда давления возрастает. При 120 = 100 мкс ведёт себя как очень пологая немонотонная функция. Все кривые для медленной волны носят убывающий характер. Очень интересно при этом отметить, что даже в случае упругого скелета (и при отсутствии межфазных взаимодействий) амплитуда давления жидкости в «медленной» волне убывает.

4. Заключение

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

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

ЛИТЕРАТУРА

1. Нигматулин Р.И. Динамика многофазных сред. Ч. 1. М.: Наука, 1987.

2. Губайдуллин А.А., Урманчеев С.Ф. Исследование прохождения волн сжатия из жидкости или газа в насыщенную пористую среду и отражения их от преград // Акустика неоднородных сред, - Вып.105, СОРАН, 1992, С. 122-128

3. Губайдуллин А.А., Дудко Д.Н., Урманчеев С.Ф. Моделирование взаимодействия воздушной ударной волны с пористым экраном // Физика горения и взрыва. 2000. Т. 36, №4. С. 87-96.

Вестник Башкирского университета.2005.№4. 19

4. Ивандаев А.И. Об одном введения «псевдовязкости» и его применения к уточнению разностных уравнений газодинамики // Журн. вычисл. математики и мат. физики. 1975. Т. 15, № 2. С. 523-527.

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

Поступила в редакцию 12.07.05 г.