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

личина давления у стенки акустической камеры примерно вдвое выше, чем в центре (эффект наложения падающей и отраженной волн). Расчетное время реверберации на частоте 1 ООО Гц составляет около 4,1с.

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

Библиографический список

1. Felippa, С. A. Partitioned analysis of coupled system/ С. A. Felippa, К. C. Park // Chapter 3 of Computational Methods for Transient Analysis, ed. by T. Belytschko and T. J. R.Hughes. North-Holland, Amsterdam, 1983. P. 157-219.

2. Орлов, А. С. Анализ нагружения сотовой панели акустическими воздействиями / А. С. Орлов // Эффек-

тивность сотовых конструкций в изделиях авиационно-космической техники : материалы II Между на р. науч,-практ. конф. Днепропетровск, 2007. С. 142-150.

3. Бате, К. Численные методы анализа и метод конечных элементов / К. Бате, Е. Вилсон. М.: Стройиздат, 1982.

4. Шимкович, Д. Г. Расчет конструкций в MSC/ NASTRAN for Windows / Д. Г. Шимкович. М.: ДМК-Пресс, 2001.

5. Лепендин, JI. Ф. Акустика / JI. Ф. Лепендин. М. : Высш. ж., 1978.

6. Виноградова, М. Б. Теория волн / М. Б. Виноградова, О. В. Руденко, А. П. Сухоруков. М.: Наука, 1979.

7. Справочник по технической акустике: пер. с нем. / под ред. М. ХеклаиХ. А. Мюллера. JI.: Судостроение, 1980.

8. Лопашев, Д. 3. Методы измерения и нормирование шумовых характеристик / Д. 3. Лопашев, Г. Л. Осипов, Е. Н. Федосеева. М.: Изд-во стандартов, 1983.

A. S. Orlov

FINITE ELEMENT MODELING OF DIFFUSIVE FIELD IN BIG VOLUME ACOUSTIC REVERBERANT CAMERAS

In this paper it is considered creation offinite element model of diffusive field for big volume acoustic reverberant cameras. Model created by usage 3d acoustic elements and different types of damping.

Keywords: modeling, diffusive field, acoustic calculation, reverberation chambers.

УДК539.3

H. В. Молокова

ЧИСЛЕННОЕ ИССЛЕДОВАНИЕ И АНАЛИЗ ПРОЦЕССА ЗАГРЯЗНЕНИЯ ПОРИСТОЙ СРЕДЫ

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

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

тая среда, углеводородный загрязнитель.

Геофильтрационная задача о движении углеводородного загрязнителя и воздуха с учетом гравитационного влияния в трехмерной постановке описывается системой уравнений:

ds

т-------Ь

dt

д(р и х) + д(р иу) + d{puz)

дх

ду

dz

= 0,

5(1-5)

т—-------- +

dt

д(р и х) + д(риу) + d(puz)

дх

ду

dz

(1)

= 0.

^ _____________*■

Щ =~ k— (V^-p^), (2, а)

h

Is f сЛ _,

и2 = ~к——(Ур2~ P2g)- (2,6)

М-2

Начальные условия:

-при/= 0: 5j=52 ,s = s0(x,y),p =р0 наГр

s = 0,p=p наГ,. (3)

’Г Г атм 2 v 7

- начальное давлениерх = рxgH распределено по

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

Граничные условия:

- на верхней границе рассматриваемой области граничное условие для s(x, у, z, I) и р(х. у, z, I) при І є (0; 7] имеет вид

5 = s0(x, у, t), р = р0 (х, у, t) на Гр (4)

р=р

L •» ЯП

Ё1.о,

дп

др др _ ф _ Ф _ Ф

dz Z = 1 9У дх У=1

k^l(-vPl -plg)= О,

K(s) =

k2(s) =

, < S < 1,

О, 0<5<5,,

3,5

(1 + 35), 0< 5< 5*,

dt dz (Xj ^ dz

+ д_кКШ ^рЛ + _^_кК^){ Фі

дх

дх

т -

5(1 -s) д

dt

дх

dz

\

АС?) Ф2

\х2 дх

АО?) Ф2

\i2 dz

-P2g

_д_

ду

АО?) Ф2

Ц2 ду

div

1K(s) | k2(s)^ h M"2 у k2(s)

^Pi-Pigk

kjjs)

h

+ div

(ft(*)-p2«)

M"2

= 0.

ницаемость ^(.v) обращается в нуль. Учитывая (7) и (8), получим

ds

т — = div

dt

АО?)

h

(va

■Рія)

(її)

(5)

где п - вектор нормали к границе Г2.

На нижней границе рассматриваемой области граничное условие для s(x,y, z, I) ир(х,у, /. I) 1(е (0; 7] имеет вид

= 0, (6)

(7)

(8)

0, 5* < 5< 1.

Математическая модель (1)...(8), основанная на двухфазной фильтрации в физических переменных с учетом гравитационного влияния, позволяет определить:

- распределение скоростей фильтрации в грунте в любой момент времени;

- распределение давления в каждой из фаз рр

- насыщенность зу

Подставим уравнения (2) в уравнения неразрывности (1). В результате получим:

дя д к^)( др1 А

т— = —к-!-------1 —-р^ +

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

Метод конечных разностей включает следующие основные этапы:

- построение сетки, охватывающей рассматриваемую область;

- построение на полученной сетке конечно-разно-стной аппроксимации, эквивалентной исходному дифференциальному уравнению и дополнительным условиям;

- формирование на основе конечно-разностной аппроксимации системы алгебраических уравнений и ее решение.

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

^,х={г,1 = /А(° ^ ^п)^а=

= гА (° - г'2 - п)>ув=

= 13Ьу(0 < /3 < п),

/,=;ф,0 <; < п), где п = ЬМ , п = ЬЛ , п = Ь,11г ,п = 77 х.

I 1 х 2 лг у 3 уР I

Заменим производные во внутренних узлах т конечно-разностными отношениями с учетом представления переменного коэффициента и получим явную разностную схему [7]:

Г.У+1 -

/1,/2,/3 /1,/2,/3

(9)

hz

Xі

/1+1/2,/2,/З

Р /1+1,/2,/3 Р /1,/2,/3

hz

-Pig

Ч І Pil,i2,B Рі1-1,і2,ІЗ п п /і-і/2,/2,/з ; Рі 8

hz

Рс(^ = Р2-Ру Уравнения (9) после преобразований примут вид (в безразмерной форме)

hx

Xі

/2+1/2,/1,/3

Рі2+\,і\,іІ Pi2,i\,il ^ hx

XJ

/2-1/ 2,/1,/3

Рі2,П,В Рі2- 1,/1,/3 ^ hx

(10)

hy

Рассмотрим важное обстоятельство, когда предельное значение насыщенности 5* = 1, при котором фазовая про-

Xі

/3+1/ 2,/1,/2

Р/3+1,/1,/2 Р/3,/1,/2

hy

Xі

/3-1/2,/1,/2

Р /3,/1,/2 Р /3—1,/1,/2

hy

SWJ§,= k ——, Swfif = k

M"2

4a:=fl* j = 0>il=0>-.,ns-l: in =0, ...,nx -1; =0, ...,«v -1;

4hh =l- h =0: 0 < i2, =3;

4,& =0> h =0; 3</2»f3 ;/V 1:

'fi i =°- tehA =3-âL =Atm- 'i =0; ÎS4À -!•

â-2(5)

si+1 \J 7

ZiL___'2L = -_JJ vi

x h

■Pi f*

ленно (в данном случае скорость продвижения фронта приблизительно равна 4 • 10'4 в безразмерных единицах), что согласуется с экспериментальными данными [9; 10]. Форму фронта определяют кривые А'(з), для этого варианта расчета в качестве модельных фазовых проницаемостей были взяты функциональные зависимости (7) и (8), а также члены дисперсионного типа. Фазовые проницаемости являются экспериментально измеряемыми функциями насыщенности. Расчеты показывают, что вблизи источника нефтезагрязнения образуется зона загрязнения, далее фронт нефтенасыщенности образуется и продвигается в основном за счет действия гравитационных сил (рис. 2).

Н

где ’kiim2(s) =X(sf + .у, , ) / 2. Разностная схема (12) аппроксимирует исходное уравнение с погрешностью 0(х + А2).

Систему уравнений (12) решаем итерационным методом. Формула позволяет получить решение s яр в явном виде. При этом на каждом временном слое вначале находим давление, а далее решаем параболическое уравнение для нефтенасыщенности. Это приводит к необходимости подбора временного шага для устойчивости счета и сходимости итерационных процессов. В численных расчетах использовались шаги сетки h = ОД, х = 0,002.

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

- пористость грунта т = 0,05.. .0,4;

- проницаемость к= 1,15 10 11,5 мкм2;

- плотность загрязнителя с 1.1... 1.4 г/см3;

- вязкость загрязнителя м = 0,54...0,9 сПз;

- высота области фильтрации загрязнителя L = 1 м;

- предельные значения насыщенности s* = 0,9; щ = 0,1;

- насыщенность на верхней границе области SQ = 1,0;

- ускорение свободного падения g = 9,8 м/с2.

Разработанная математическая модель фильтрации

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

Для проведения численных расчетов была реализована программа с использованием технологии визуального программирования Borland Delphi [8].

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

Рис. 1. Динамика распространения фронта при различных значениях безразмерного времени t 408; 857; 1 308

Рис. 2. Распределение давления при различных значениях безразмерного времени Г: 408; 857; 1 308

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

Полезные качественные сведения о закономерностях движения углеводородов в пористой среде могут быть получены путем анализа природных условий и постановки экспериментов. Проведенные лабораторные исследования пропитки воздушно-сухого песка (II' 10... 12 %) при комнатной температуре различными объемами нефти (Г ’ = 150 мл и |’ = 100 мл) показали, что процесс протекания нефти происходит преимущественно в вертикальном направлении, но наблюдается также некото-

рое растекание в боковые стороны. Эти данные совпадают с экспериментальными результатами, описанными в [9].

*

Рис. 3. Изолинии нефтенасьпценности Г = 1 500

Было установлено, что вся нефть протекла в хорошо проницаемый песок со скоростью, приблизительно равной

3,7 • Ю ' м/ч. Это объясняется тем, что гравитационный напор оказывается сильнее капиллярного противодействия в течение всего времени наблюдения. Чем меньше высота нефтяного слоя, тем слабее его растекание в стороны, поэтому боковые потоки ослабевают книзу (рис. 4).

Рис. 4. Динамика проникновения нефти при начальной высоте нефтяного слоя Ь( = 10 мм: Н - высота области фильтрации; Ь - полуширина области фильтрации;

В = 113 мм - ширина углеводородного пятна на поверхности; V = 100 мл - обьем нефти

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

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

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

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

Библиографический список

1. Численные методы задач математической физики / Н. С. Бахвалов. Г. М. Кобельков, Ю. А. Кузнецов и др. // Современные проблемы вычислительной математики и математического моделирования. М. : Наука, 2005. Т. 1. С. 18-28.

2. Костерин, А. В. Моделирование загрязнения почв и грунтов органическими жидкостями / А. В. Костерин // Современные проблемы математического моделирования : тр. Всерос. школы-семинара. Ростов н/Д, 2001. С. 142-147.

3. Булгакова, Г. Т. Неравновесная двухфазная фильтрация / Г. Т. Булгакова. А. В. Жибер, Т. А. Файдуллин // Математическое моделирование РАН. 2006. Т. 18. № 10. С. 19-38.

4. Самарский, А. А. Математическое моделирование. Процессы в сложных экономических и экологических системах/ А. А. Самарский, Н. Н. Моисеев, А. А. Петров М.: Наука, 1986.

5. Гринин, А. С. Математическое моделирование в экологии : учеб. пособие для вузов / А. С. Гринин. М. : ЮНИТИ-Дана, 2003.

6. Ляшко, И. И. Численное решение задач тепло и мас-сопереноса в пористых средах / И. И. Ляшко, Л. И. Демченко. Киев: Наук, думка, 1991.

7. Тихонов, А. Н. Разностные методы решения многомерных задач / А. Н. Тихонов, А. А. Самарский. М.: Наука, 1989.

8. Бежанова, М. М. Практическое программирование. Визуальное программирование в среде □е1р1п : учебник /М. М. Бежанова, Л. А. Москвина. М.: Логос, 2001.

9. Грищенко, А. И. Экология. Нефть и газ / А. И. Гри- 10. Басниев, К. С. Подземная гидромеханика / К. С. Бас-

щенко, Г. С. Акопова, В. М. Максимов; ВНИИЭГазпром. ниев, И. Н. Кочина, В. М. Максимов. М.: Недра, 1993.

М., 1995.

N. V. Molokova

COMPUTATIONAL INVESTIGATION AND THE PROCESS OF THE POROUS MEDIUM CONTAMINATION ANALYSIS

It is covered a two-dimensional model of filtering, taking into account the movement ofair and hydrocarbon pollutants into the porous soil. The model includes a system of partial differential equations with additional conditions. Among the differential equations balance equation element in a porous medium - inseparability equation, and differential equations of motion are included. For circuit system the equations of the considered pollutant state and environment. The are entered initial and boundary conditions correspond to the filtration process, beginning with the ground surface and initial stage of oil products ladling. A comparative analysis of the results of mathematical modeling with experiments is carries out.

Keywords: two-dimensional model of filtering, geofiltration, mathematical modeling, porous medium, hydrocarbon pollutants

УДК519.688

К. В. Богданов, А. Н. Ловчиков

МОДЕЛИРОВАНИЕ СИСТЕМ С СУЩЕСТВЕННЫМИ НЕЛИНЕЙНОСТЯМИ НА ОСНОВЕ ПАРАЛЛЕЛЬНЫХ ВЫЧИСЛЕНИЙ С АСИНХРОННЫМ ОБМЕНОМ

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

Ключевые слова: моделирование, параллельные вычисления, EDA-системы.

Современные EDA-системы развивались в течение длительного времени. Но это развитие было экстенсивным и во многом консервативным, это обусловлено следующими причинами: во-первых, необходимость сохранения обратной совместимости, т. к. объем накопленной информации для EDA- и CAE-систем весьма велик, и по сей день используются различные версии систем; во-вто-рых, на момент зарождения наиболее известных систем, таких как MATLab, Spice и т. п., средства разработки и существующее аппаратное обеспечение диктовали классический подход к проектированию программных систем, при котором имеется единственный вычислительный поток, в котором и производится расчет по заданной математической модели.

Рассмотрим одну из типичных задач, ставящихся перед специализированными ED A-системами, проведение анализа переходных процессов для электронной схемы. Обычная EDA-система, такая как, к примеру, Microcap, организует этот процесс в несколько этапов (рис. 1).

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

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

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

Основной этап - решение системы обыкновенных линейных дифференциальных уравнений (ОДУ) численным методом. Начальные условия берутся из параметров, задаваемых пользователем перед проведением процедуры анализа. На этом этапе для упрощения вычислений зачастую проводится процедура линеаризации.

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