Бадалова М.В. Рефлексивная организация этического пространства п. Рефлексивная организация этического пространства психолога
 Скачать 1.31 Mb.
|
|
Прозорова Гульшат Ринатовна старший преподаватель Сургутского государственного педагогического университета, г. Сургут OF THE NUMERICAL SOLUTION OF PROBLEMS OF THE TWO- PHASE FILTRATION ON THE RELATED MOBILE GRIDS Prozorova Gulshat assistant of Surgut State Pedagogical University, Surgut АННОТАЦИЯ В данной статье рассматривается решение двумерных задач двухфазной фильтрации на подвижных сетках, пространственная структура которых определяется самими значениями рассчитываемого решения. Для задач нелинейной теплопроводности применение таких сеток, представляющих собой совокупность подвижных изотерм, дает эффективный способ выделения и прослеживания таких особенностей, как движущиеся границы фазовых переходов и бегущие температурные волны. В статье дан алгоритм численного метода к решению двумерных задач фильтрации на подвижных сетках, привязанных к значениям рассчитываемого поля температур, с обобщением на класс задач многофазной фильтрации paper deals with the two-dimensional problems of two-phase filtration movable grids, the spatial structure is determined by the values calculated decision. For the purposes of the application of the nonlinear heat grids, representing an aggregate of moving isotherms, provides an effective method for isolating and tracking features such as moving boundaries of phase transitions and running temperature waves. In the article the algorithm for the numerical method to solve the two-dimensional filtration problems on moving grids linked to the values calculated temperature field, with the generalization of the class of problems of multiphase Ключевые слова численное моделирование, двухфазная фильтрация, слоисто–неоднородные пористые среды, вытеснение нефти, пространственная сетка numerical modeling, two-phase filtration, layered porous medium, oil displacement, spatial Для исследования характера многофазных фильтрационных течений в неоднородных пластах широко применяются методы математического моделирования. Решение практических задач многофазной фильтрации возможно только с применением численных методов. Проблемы построения эффективных вычислительных алгоритмов для данного достаточно сложного класса нелинейных задач хорошо известны и далеко не преодолены [2], [3], Рассмотрим совместное течение двух несжимаемых жидкостей в среде с постоянной пористостью и проницаемостью. Полагая в качестве неизвестных s насыщенность порового пространства первой фазой и P – давление во второй фазе, система двух уравнений относительно определяемых полей насыщенности и давления в безразмерном виде имеет следующий вид ) ( ) ( ) ( ) ( ) σ σ σ σ ∂τ ∂σ grad div grad div 1 1 J k Ï P k ′ − = (1) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 0 grad div grad div 1 2 1 = ′ − ⋅ + σ σ σ σ µ σ J k Ï P k k (2) 54 Евразийский Союз Ученых (ЕСУ) # 2 (23), 2016 | ФИЗИКО - МАТЕМАТИЧЕСКИЕ НАУКИ где P P P P ∆ − = 0 – приведенное нормированное давление (P0 и DP – характерные значения давления и его перепада, 2 / l at = τ – безразмерное время (l – пространственный масштаб отношение вязкостей двух фаз, k m P Ï n 0 cos θ α = – безразмерный параметр, характеризующий влияние капиллярных сил. Система (1), (2) дополняется стандартными граничными условиями для давления и граничными и начальным условиями для насыщенности, определяющимися конкретной постановкой задачи. Подходы, лежащие в основе предлагаемого численного алгоритма, рассмотрим на примере решения уравнения вида ) grad ) ( ( div ) ( σ σ τ σ σ k a = ∂ ∂ (3) содержащего все дифференциальные операторы, входящие в систему (1), (Пусть уравнение (3) решается в области } ), , ( ) , ( : ) , ( { 0 0 K N z z D η η η τ η ξ τ η η ξ ≤ ≤ ≤ ≤ = в некоторой системе координат (x, h). Допускается, что границы областей x = z0(h,t) и x = zN(h,t) могут перемещаться стечение времени. Для дискретизации расчетной области используется полукриволинейная подвижная сетка } ..., ,1 , 0 ; ..., ,1 , 0 : ) , ( { K k N n D k nk = = = η ξ , составленная совокупностью узлов, соответствующих точкам пересечения семейства фиксированных координатных линий h = hk и семейства подвижных линий вида x = zn(h,t). Для последних должно выполняться требование ) , ( ) , ( ) , ( 1 1 τ η τ η τ η + − ≤ ≤ n n n z z z для любого n и для каждой пары значений (h,t). Фрагмент такой расчетной сетки представлен на рис. На каждый момент времени с любым узлом (n, k) связана тройка величин (xnk(t), hk, snk(t)), где первая из координат) и значение решения в этом узле snk(t) = s(xn(t), hk) могут изменяться во времени. Положение узла на координатной линии h = hk может определяться заданием одной из величин координаты xnk(t) или значения узловой величины snk(t). Таким образом, расчетная сетка может включать узлы двух типов перемещающиеся вдоль данной координатной линии по заданному закону и прослеживающие перемещение точки с заданным изменением рассчитываемой узловой характеристики. В частном случае величины и snk могут быть фиксированными. Рисунок 1. Фрагмент расчетной сетки В рамках сделанных допущений можно определить расчетные сетки самой разнообразной структуры, сочетающей узлы двух описанных типов. Построение системы уравнений, представляющих дискретный аналог исходного дифференциального уравнения (3), основывается на интегральных тождествах, справедливых для окрестности каждого узла расчетной сетки. Пусть (n, k) – произвольный внутренний узел. Его сеточная окрестность D(n, k) представляется объединением элементарных подобластей 1 , 0 , , ) , ( = − − = l j l k j n k n D D где } ), , ( ) , ( : ) , ( { 1 1 + + ≤ ≤ ≤ ≤ = k k n n nk z z D η η η τ η ξ τ η η ξ (рис. Интегральное тождество для узла получается интегрированием уравнения (3) по каждой из подобластей, составляющих его окрестность, со своими весовыми функциями и последующим суммированием результатов интегрирования. Веса интегрирования определяются с помощью набора функций следующего вида Евразийский Союз Ученых (ЕСУ) # 2 (23), 2016 | ФИЗИКО - МАТЕМАТИЧЕСКИЕ НАУКИ 55 , , , 1 1 ) ( 1 1 ) ( n n n n n n n n n z z z z z z z − = ∆ ∆ − = ∆ − = + + + − − − ξ ν ξ ν , , , 1 1 ) ( 1 где нижний индекс совпадает с индексом узла, а верхний определяет ориентацию подобласти относительно данного узла. Для подобласти Dnk узла (n, k) весовая функция есть произведение ) ( ) ( + + k n v ω , для подобласти Dn-1,k этого же узла – произведение ) ( ) ( + − k n v ω и т.д. Далее будем считать, что координатная система (x, h) является ортогональной, так что развернутое представление уравнения (3) имеет вид ∂ ∂ ∂ ∂ + ∂ ∂ ∂ ∂ = ∂ ∂ η σ σ η ξ γ η ξ σ σ η ξ β ξ σ σ η ξ α ) ( ) , ( ) ( ) , ( ) ( ) , ( k k t a (4) где коэффициенты a(x, h), b(x, h) и g(x, h) – известные определяющиеся координатной системой функции Результат интегрирования уравнения (4) по подобласти Dnk узла (n, k) приводится к равенству где η ξ ν τ σ ∂ σ η ξ α ω ξ ξ η η d d a I n k k n n n k k ) ( ) ( 1 1 ) ( ) , ( + + ∫ ∫ + + ∂ = η ξ ∂ξ ∂σ σ η ξ β ω η η d d k z F n n k k z z n k k n ∫ ∫ + + ∆ = + 1 1 ) ( ) , ( ) ( η ξ ∂η ∂ν ∂η ∂σ σ η ξ γ ω η η d d k Y n z z k k n n n k k ) ( ) ( 1 1 ) ( ) , ( + + ∫ ∫ + + = η ξ ν ∂η ∂σ σ η ξ γ η η η d d k K n z z k k n n n k k ) ( 1 1 ) ( ) , ( 1 + ∫ ∫ + + ∆ = 0 ) , ( ) , ( ) ( ± = ± − − = n z n z z k Q ξ ∂η ∂ ∂η ∂σ η ξ γ ∂ξ ∂σ η ξ β σ ) ( ) , ( 0 ± = ± − = k k Q η η η ∂η ∂σ σ η ξ γ , Для остальных подобластей окрестности узла получаются равенства, аналогичные (5). Заметим, что правые части равенств (5) представлены потоковыми членами по общим границам подобластей и при суммировании взаимно уничтожаются. В итоге интегральное тождество по окрестности произвольного внутреннего узла будет иметь вид ( ) 0 ) 1 ( ) 1 ( 1 , 0 , , 1 = − + + − + ∑ = − − + l j l k j n l j K Y F I (6) Для узлов, лежащих на границе расчетной области, тождества интегрального баланса получается аналогично за исключением того, что окрестность граничного узла состоит из одной или двух подобластей, а члены, выражающие потоки на внешних границах области заменяются заданными функциями из граничных условий. Система численных уравнений относительно неизвестных узловых величин выводится на основе тождеств (6) последовательным раскрытием интегралов по обеим координатам. Для неизвестного распределения s(x, h) в пределах каждой подобласти используется линейная аппроксимация по координате x вдоль каждого фиксированного значения координаты где ) , ), , ( ( ) , ( τ η τ η σ τ η σ n n z = – распределение поля s вдоль линии Для вычисления интегралов по координате h с весами wk(h) применяются квадратурные формулы 1 )) ( 2 ) ( ( 6 ) ( )), ( 2 ) ( ( 6 ) ( 1 ) ( 1 являющиеся следствием линейной аппроксимации подынтегральной функции f(h) в пределах отрезков [hk, При получении аппроксимирующих уравнений входящие в (4) коэффициенты в пределах каждой элементарной подобласти Dnk заменяются их осредненными по подобла- сти значениями. Производные от узловых величин повремени аппроксимируются неявными разностными отношениями, гарантирующими безусловную устойчивость вычислительной схемы от величины временного шага, что является принципиально важным при использовании расчетных сеток с большой пространственной неравномерностью. В конечном итоге расчет одного временного шага приводит к девятиточечной системе алгебраических уравнений относительно неизвестных величин в точках расчетной сетки, входящих в единую узловую окрестность. В случае явного задания расположения узлов xnk получаемая система уравнений линейна. В общем случае, когда сетка включает узлы с заданными значениями определяемого решения snk, система является нелинейной относительно координат xnk. Для ее решения используется линеаризация по методу Ньютона, приводящая на каждой итерации к линейным девятиточечным уравнениям с блочно-трехдиагональной матрицей. К последним применяется матричная прогонка [5]. Отметим, что численная реализация алгоритма не зависит от конкретной структуры расчетной сетки. При этом могут использоваться различные координатные системы, а также разбиение расчетной области наконечное число подобластей, геометрия которых для подходящей системы координат не противоречит оговоренным соглашениям. В качестве примера рассматривается процесс вытеснения нефти водой в однородном пласте на элементе пяти- точечной системы скважин. С учетом симметрии задачи исследуемая область представляет единичный квадрат с круговыми вырезами радиуса r0 с центрами в двух противоположных углах, соответствующих контурам нагнетающей и добывающей скважин: Постановка задачи включает уравнения двухфазной фильтраци с граничными условиями: на контуре нагнетающей скважины 0 , 0 , 2 0 2 2 > > = + y x r y x : 1 , 1 = = σ P на контуре добывающей скважины 1 ,1 , ) 1 ( ) 1 ( 2 0 2 2 < < = − + − y x r y x : 0 = P ; ) , ( ) , ( , ) , ( ) , ( ) , , ( 1 ) ( 1 1 ) ( τ η ξ τ η τ η σ τ η σ τ η ξ σ + − + + + ≤ ≤ + ≈ n n n n n n z z v v } 1 0 , 1 0 , )1 ( )1 ( , :) , {( 2 0 2 2 2 0 2 2 ≤ ≤ ≤ ≤ ≥ − + − ≥ + y x r y x r y x y x 56 Евразийский Союз Ученых (ЕСУ) # 2 (23), 2016 | ФИЗИКО - МАТЕМАТИЧЕСКИЕ НАУКИ условиями непротекания на остальных границах области начальным условием Для численного решения использовалась координатная система ϕ η ρ ξ = = , ln , где r и j – радиус и угол полярных координат. Расчет поля давлений и поля насыщенности производился на отдельных сетках. Для учета больших градиентов давления в прискважин- ных зонах и получения пространственной дискретизации, близкой к конфигурации изобар и линий тока удобнее разбить исследуемую область на две подобласти D1 и D2 диагональю проходящей через вершины, в которых нет скважин, и перейти для каждой подобласти к своим полярным координатам. В каждой из двух подобластей полярный полюс располагаем в той вершине, где находится ось скважины. Таким образом, } 2 / 0 , ) ( ln ln : ) , {( 0 2 1 π η ϕ ξ η ξ ≤ ≤ ≤ ≤ = = g r D D , где ) sin /(cos 1 ) ( ϕ ϕ ϕ + = g . В каждой области задавалась равномерная по координатами расчетная сетка, дающая автоматическое сгущение узлов в узких прискважинных зонах, где поле давлений имеет большие градиенты Поле насыщенности рассчитывалось в системе координат) с полюсом в центре нагнетающей скважины на подвижной сетке, соответствующей совокупности линий x = zn(h,t) с одинаковыми значениями насыщенности sn вин- тервале snÎ[0; 1]. Таким образом, расчет распределения s(x, h, t) сводился к прослеживанию перемещения вдоль прямых h = hk точек линий изонасыщенности с заданными значениями. Для интервала изменения величин s и h задавалось равномерное разбиение. Положение изолинии x = z0(h,t) с насыщенностью s0=1 фиксировано и совпадает с контуром скважины. Изолиния x = zN(h,t) с насыщенностью sN=0 представляет подвижную границу области двухфазного течения с соответствующим граничным условием Связанная система уравнений (5), (6) является внешне нелинейной относительно определяемых полей давления и насыщенности. Такой же нелинейностью характеризуется и система численных уравнений при использовании полностью неявной аппроксимации повремени. Поэтому расчет каждого шага повремени требует внешних итераций. Последовательность действий для перехода на новый временной слой включала следующие этапы) Для заданного распределения насыщенности из дискретного аналога уравнения (2) находится поле давлений. Уравнение (2) является линейным относительно давлений и решается на фиксированной сетке, поэтому система численных уравнений также линейна) Исходя из найденного поля давлений из дискретного аналога уравнения (1) определяется поле насыщенности. Так как для нелинейного уравнения (1) используется подвижная сетка с фиксированными значениями насыщенности, то система численных уравнений нелинейна относительно координат узлов и решается итерациями по методу Ньютона Оба этапа итерируются до достижения сходимости внешних итераций по насыщенности. Так как каждое из определяемых полей рассчитывается на разных пространственных сетках, то перед каждым этапом производится интерполяция давления и насыщенности с одной сетки на другую. Предлагаемый численный алгоритм был протестирован на рассмотренной задаче с использованием последовательности сгущающихся сеток для различных значений безразмерных величин µ и Ï . Тестирование показало, что хорошее пространственное разрешение особенностей рассчитываемых полей достигается на сетках с небольшим числом узлов. В качестве примера на рис. 2 представлены поля давления и насыщенности, рассчитанные с помощью изложенного алгоритма на сетках 16´16 для µ = 0.1 и Ï = Рисунок 2. a) Распределение давления b) Распределение насыщенности Численное моделирование задачи двухфазной фильтрации позволяет использовать созданный в ней математический, алгоритмический и программный продукт для решения задач прогнозирования многофазных фильтрационных течений, в частности в пластах сложной структуры. Результаты численного исследования дают ряд новых выводов об особенностях вытеснения нефти из слоисто-неоднородной среды, которые могут найти практическое применение при совершенствовании технологий разработки нефтяных пластов с послойной неоднородностью. Список литературы. Влияние межпластовых перетоков и капиллярных сил на процесс вытеснения нефти в слоисто-неоднородном пласте Текст / Ю. А. Сигунов, ГР. Усманова // Известия РАН. 0 0 ) ( = ∂ ∂ ∂ ∂ − ∂ ∂ + = η ξ η η σ ξ σ N z N z Евразийский Союз Ученых (ЕСУ) # 2 (23), 2016 | ФИЗИКО - МАТЕМАТИЧЕСКИЕ НАУКИ Механика жидкости и газа. - 2007. - N 6. - С. 85-92 : 4 рис. - Библиогр.: с. 92 (7 назв. ) 2. Математическое моделирование пластовых систем : перс англ. / Х. Азиз, Э. Сеттари . – е изд, стереотип . – М. : Ин-т компьют. исслед., 2004 . – 416 с. – (Современные нефтегазовые технологии) . - ISBN 5-939723-55-1 3. Коновалов АН. Задачи фильтрации многофазной несжимаемой жидкости. – Новосибирск Наука, 1988. — 166 с. — ISBN 5-02-028569 -2. 4. Сигунов Ю.А. Методы решения классической задачи Стефана – Сургут РИО Сургутского государственного педагогического университета, 2009. – с. 5. Самарский А.А., Гулин А.В. Численные методы Учебное пособие для вузов. – М Наука. Гл. ред. физмат. лит, 1989. – 432 с. Флетчер К. Вычислительные методы в динамике жидкостей. В х томах. – М Мирт (т TWO GRAPH TRANSFORMATION |