3.2 Фильтрация однофазной несжимаемой жидкости

3.2.1 Постановка

Рассмотрим фильтрацию несжимаемой жидкости в трехмерной области Ω:

Ω = {(x,y,z), lmxin < x < lmxax, lmyin < y < lmyax, lmzin < z < lmzax}.
(3.3)

Течение жидкости в пористой среде описывается уравнением сохранения массы и законом Дарси, для удобства отсчет координаты z возьмем вниз по вертикали:

- ∇ (ρu ) = ∂-(ρϕ)+ ρq,
         ∂t
(3.4)

     k-
u = -μ (∇p + ρg)
(3.5)

где ρ – плотность жидкости, u – скорость течения, ϕ – пористость области, q – функция источников закачки и откачки, k – тензор проницаемости среды, μ – вязкость жидкости, g – вектор ускорения свободного падения.

Предположим, что вязкость μ, проницаемость k и пористость ϕ не зависят от давления и пространственных переменных. Поскольку рассматриваем несжимаемую жидкость, то плотность ρ = const. Подставляем уравнение скорости (3.5) в уравнение сохранения массы (3.4) и делим на плотность ρ. Тогда получаем следующее линейное уравнение для давления в области Ω:

  (      )
∇  - k∇p   = q, (x,y,z) ∈ Ω.
     μ
(3.6)

Левая часть полученного уравнения представляет собой оператор диффузии с постоянным коэффициентом -k
μ. Дополним граничными условиями первого рода по горизонтали, т.е. по координатам x и y, и условием не протекания по вертикали:

p(lmxin,y,z) = p(lmxax,y,z) = p0 + ρg(z - lmzin),
p(x,lmyin,z) = p(x,lmyax,z) = p0 + ρg(z - lmzin),
(3.7)

  |          |
∂p||      = ∂p||      = - ρg.
∂z|z=lmzin   ∂z|z=lmzax
(3.8)

3.2.2 Аппроксимация

Для численного решения поставленной задачи определим равномерную прямоугольную сетку:

    (       |     min          --------          max    min  )
    {       ||x = lx  + ihx, i = 1,Nx---1, Nxhx = lx   - lx  , }
ω = ((x,y,z)||y = lmyin + jhy, j = 1,Ny---1,  Nyhy = lmyax - lmyin, ) .
            |z = lzmin + khz, k = 1,Nz - 1, Nzhz = lmzax - lmzin
(3.9)

Тогда уравнение (3.6) аппроксимируем для внутренних узлов следующим образом:

d  1  pi+1,j,k-2-pi,j,k - d  1  pi,j,k-p2i-1,j,k+
 i+ 2,j,kpi,j+h1,xk-pi,j,k    i- 2,j,kpi,j,kh-xpi,j-1,k
di,j+ 12,k    h2y     - di,j- 12,k    h2y     +      , (x,y,z) ∈ ω.
di,j,k+12 pi,j,k+h1-2z-pi,j,k - di,j,k- 12 pi,j,k-hp2zi,j,k-1 = qi,j,k
(3.10)

где di,j,k – коэффициенты для разностного оператора диффузии. Граничные условия (3.7) и (3.8) аппроксимируются в следующие конечные разности:

p0,j,k = pNx-1,j,k =
pi,0,k = pi,Ny-1,k = p0 + ρg(khz - lmzin),
 pi,j,1-hpzi,j,0-= - ρg, pi,j,Nz-hpzi,j,Nz-1-= - ρg.
(3.11)

Подставляем аппроксимацию граничных условий (3.11) в систему уравнений (3.10), тогда коэффициент di,j,k в полуузлах по координатам x и y принимает следующие значения:

d  1   = - k, i = 1,N---1,
 i2,j,k    μk     ---x----
di,j12,k = - μ, j = 1,Ny - 1.
(3.12)

А для граничных полуузлов по z учитывается условие Неймана (3.8):

d     1= 0, k = 1,
di,j,k- 21= - k, k = 2,N---2,
 i,j,k 2    μ        z
di,j,k+ 12 = 0, k = Nz - 1.
(3.13)

Оставшаяся часть переходит в правую сторону и получаем следующую систему уравнений:

Dp = f, (x,y,z) ∈ ω,
(3.14)

где D – разностный оператор диффузии с коэффициентами (3.123.13). Правая часть f без учета граничных узлов принимает следующие значения:

fi,j,k = qi,j,k.
(3.15)

На граничных узлах, при i = 1 в правую часть f добавляется выражение:

- d1  (p0 + ρg(khz - zmin))
---2,j,k------2------------,
           hx
(3.16)

а при i = Nx - 1,j = 1,j = Ny - 1,k = 0 и k = Nz - 1 добавляются соответственно:

- d   1  (p0+ρg(khz-zmin))
--Nx--2,j,k--h2x---------- ,
  -di,12,k(p0+ρg(khz-zmin))-
           h2y         ,
--di,Ny- 12,k(p0+2ρg(khz-zmin)),
       - d hy1(ρg)
       --i,j,hz2-- ,
      -di,j,Nz- 12(ρg)
           hz     .
(3.17)

3.2.3 Реализация

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

// Parameters.h 
class Parameters { 
public: 
    Parameters(); 
 
    double permeability; 
    double muFluid; 
    double rhoFluid; 
    double g; 
 
    double p0; 
}; 
    

В реализации класса Parameters зададим значения по-умолчанию в конструкторе:

// Parameters.cpp 
#include "Parameters.h" 
 
Parameters::Parameters() { 
    double darcy  = 9.869233 * 1.0e-13; 
    double mDarcy = 1.e-3 * darcy; 
    permeability  = 50 * mDarcy; 
    muFluid       = 0.001; 
    rhoFluid      = 1000.0; 
    g             = 9.8; 
 
    p0            = 1.e7; 
} 
    

Определим разностный оператор диффузии для давления унаследовав класс DiffusionOperator:

// PressureDiffOperator.h 
#include "DiffusionOperator.h" 
#include "Parameters.h" 
 
class PressureDiffOperator : public DiffusionOperator { 
public: 
    PressureDiffOperator(Parameters parameters); 
    virtual ~PressureDiffOperator(); 
 
    virtual void prepare(int level); 
 
protected: 
    virtual double diff(int i, int j, int k, int di, int dj, int dk); 
 
private: 
    Parameters p; 
 
}; 
    

В реализации класса PressureDiffOperator задаем коэффициент диффузии из формул (3.12 - 3.13):

// PressureDiffOperator.cpp 
#include "PressureDiffOperator.h" 
 
PressureDiffOperator::PressureDiffOperator(Parameters parameters) 
: DiffusionOperator() 
, p(parameters) { 
} 
 
PressureDiffOperator::~PressureDiffOperator() { 
} 
 
void PressureDiffOperator::prepare(int level) { 
    DiffusionOperator::prepare(level); 
} 
 
double PressureDiffOperator::diff(int i, int j, int k, int di, int dj, int dk) { 
    if ((k == 0          && dk == -1) || 
        (k == (mN.z - 1) && dk ==  1)) 
        return 0.0; 
 
    return - p.permeability / p.muFluid; 
} 
    

Правую часть системы уравнений (3.14) определим через унаследование класса Formula:

// PressureRhsFormula.h 
#include "Formula.h" 
#include "Variable.h" 
#include "Parameters.h" 
 
class PressureRhsFormula : public Formula { 
public: 
    PressureRhsFormula(Parameters parameters, Variable& well); 
    virtual ~PressureRhsFormula(); 
 
    virtual void prepare(int level); 
    virtual double function(int i, int j, int k); 
 
 
private: 
    string mWellName; 
    DAValue* mWell; 
 
    Parameters p; 
}; 
    

Реализация класса PressureRhsFormula задается следующим образом:

// PressureRhsFormula.cpp 
#include "PressureRhsFormula.h" 
 
PressureRhsFormula::PressureRhsFormula(Parameters parameters, Variable& well) 
: Formula() 
, p(parameters) 
, mWellName(well.getName()) { 
} 
 
PressureRhsFormula::~PressureRhsFormula() { 
} 
 
void PressureRhsFormula::prepare(int level) { 
    Formula::prepare(level); 
 
    mWell = &(Storage::instance().get(mWellName).getDAValue(level)); 
} 
 
double PressureRhsFormula::function(int i, int j, int k) { 
    double qTotal = mWell->get(i, j, k); 
    double result = qTotal / mH.x / mH.y / mH.z; //*h.x*h.y 
 
    if (i == 0) { 
        result += p.permeability / p.muFluid * (p.p0 + p.rhoFluid * p.g * mH.z * k) / mH.x / mH.x; 
    } 
    if (j == 0) { 
        result += p.permeability / p.muFluid * (p.p0 + p.rhoFluid * p.g * mH.z * k) / mH.y / mH.y; 
    } 
    if (i == mN.x - 1) { 
        result += p.permeability / p.muFluid * (p.p0 + p.rhoFluid * p.g * mH.z * k) / mH.x / mH.x; 
    } 
    if (j == mN.y - 1) { 
        result += p.permeability / p.muFluid * (p.p0 + p.rhoFluid * p.g * mH.z * k) / mH.y / mH.y; 
    } 
    if (k != (mN.z - 1)) { 
        result += p.permeability / p.muFluid * p.rhoFluid * p.g / mH.z; 
    } 
    if (k != 0) { 
        result -= p.permeability / p.muFluid * p.rhoFluid * p.g / mH.z; 
    } 
 
    return result; 
}

Теперь осталось задать формулу для скважин и основную логику программы. Описание класса скважин WellFormula:

// WellFormula.h 
#include "Formula.h" 
 
enum WellType { 
    WELL_WETTING, 
    WELL_NONWETTING 
}; 
 
class WellFormula : public Formula { 
public: 
    WellFormula(double yield, WellType type); 
    virtual ~WellFormula(); 
 
    virtual void prepare(int level); 
    virtual double function(int i, int j, int k); 
 
private: 
    double mYeild; 
    WellType mWellType; 
}; 
    

Зададим две вертикальные скважины по всей области:

// WellFormula.cpp 
#include "WellFormula.h" 
 
WellFormula::WellFormula(double q, WellType type) 
: Formula() 
, mQTotal(0.002) 
, mWellType(type) { 
} 
 
WellFormula::~WellFormula() { 
} 
 
void WellFormula::prepare(int level) { 
    Formula::prepare(level); 
} 
 
double WellFormula::function(int i, int j, int k) { 
    int wellHeight = mN.z; 
    double result = 0.0; 
 
    if ((((i == mN.x * 1 / 2) && (j == mN.y * 1 / 3)) || 
         ((i == mN.x * 1 / 2) && (j == mN.y * 2 / 3))) && 
        (k <= wellHeight)) { 
 
        if (mWellType == WELL_WETTING) { 
            result = -mQTotal / wellHeight; 
        } else { 
            result = mQTotal / wellHeight; 
        } 
    } 
 
    return result; 
} 
    

Основное тело программы:

// main.cpp 
#include "CompositeCommand.h" 
#include "FormulaSolver.h" 
#include "LinearSolver.h" 
#include "InformCommand.h" 
#include "SaveCommand.h" 
#include "Logger.h" 
#include "AdditiveEquation.h" 
 
#include "PressureDiffOperator.h" 
#include "PressureRhsFormula.h" 
#include "WellFormula.h" 
 
void simulate() { 
    Parameters params; 
 
    Variable well("well"); 
    WellFormula wellFormula(0.002, WELL_WETTING); 
 
    Variable pressure("pressure"); 
    PressureDiffOperator pDiff(params); 
    PressureRhsFormula   pRhs(params, well); 
 
    AdditiveEquation pEquation(pRhs); 
    pEquation.addOperator(pDiff); 
 
    CompositeCommand problem; 
        FormulaSolver wellSolver(well, wellFormula); 
        problem.addChild(wellSolver); 
 
        LinearSolver pressureSolver(pressure, pEquation); 
        problem.addChild(pressureSolver); 
        InformCommand pInfoMax(pressure, InformMax); 
        problem.addChild(pInfoMax); 
        InformCommand pInfoMin(pressure, InformMin); 
        problem.addChild(pInfoMin); 
        SaveCommand pSave(pressure); 
        problem.addChild(pSave); 
 
    problem.execute(); 
 
    Storage::instance().clear(); 
} 
 
int main(int argc, char** argv) { 
    PetscInitialize(&argc, &argv, (char *) 0, ""); 
    simulate(); 
    PetscFinalize(); 
    return 0; 
} 
    

3.2.4 Численный эксперимент

При расчетах использовались следующие входные параметры модели: k = 0.05 Darcy, ϕ = 0.2, μ = 0.001 Pa s. На рисунке (3.8) приводится график среза распределения давления в области lxmin = lymin = lzmin = 0, lxmax = lymax = 500 m, lzmax = 100 m. Некоторые из параметров задаются через конфигурационный файл main.config:

-log_file main.log 
-log_priority DEBUG 
-log_to_stdout true 
-outdir ./out/ 
 
-ksp_type gmres 
-pc_type jacobi 
 
-grid_levels 1 
 
-grid_nx 21 
-grid_ny 21 
-grid_nz 21 
 
-grid_minx 0.0 
-grid_miny 0.0 
-grid_minz 0.0 
 
-grid_maxx 500.0 
-grid_maxy 500.0 
-grid_maxz 100.0 
 
-value_pressure 1.0e7 
        

Конфигурационный файл задается через консольный параметр -config при запуске программы.


PIC
Рис. 3.8: Срез распределение давления при y = 250 m.









1 процесс2 процесса4 процесса8 процессов16 процессов






106 20.13 15.20 6.91 3.27 1.63






2 106 58.38 45.47 22.98 10.46 4.91






4 106 138.85 109.35 54.58 27.62 12.35






8 106 282.25 227.36 111.70 55.41 27.96







Таблица 3.1: Зависимость времени счета (сек.) от количества неизвестных и запущенных процессов.