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

3.3.1 Постановка

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

Рассмотрим основные законы, описывающие двухфазную фильтрацию. Запишем законы сохранения массы для водной и нефтяной фазы:

∂-(ϕρwsw)+ div(ρwuw) = ρwqw,
∂t
-∂(ϕρoso)+ div(ρouo) = ρoqo,
∂t
(3.18)

Закон Дарси для каждой фазы:

u α = - kkrαgrad(pα + ραgz), α = o,w.
        μα
(3.19)

Считаем, что фазы полностью заполняют поровое пространство

so + sw = 1.
(3.20)

Фазовые давления связаны через капиллярное давление

pc = po - pw.
(3.21)

В уравнении (3.18) распишем производные по времени и получим:

ϕρ ∂sw-+ s ∂-(ϕρw-)∂p-+∇ (ρ u  ) = ρ q ,
  w ∂t    w  ∂p   ∂t      w w     w w
ϕρo∂so + so∂(ϕρo)∂p-+∇ (ρouo) = ρoqo.
    ∂t       ∂p  ∂t
(3.22)

Для получения уравнения для давления, домножим каждое уравнение на ρ1
 α и сложим с учетом соотношения (3.20)

 ∑   sα ∂(ϕρα)∂p   ∑    1            ∑
     ρ----∂p--∂t-+      ρ-∇ (ραu α) =     qα.
α=o,w  α            α=o,w  α          α=o,w
или
 ∑   sα∂(ϕρα) ∂p         ∑    uα       ∑
     ρ---∂p-- ∂t + ∇uT +      ρ-∇ ρα =     qα.
α=o,w  α                  α=o,w  α      α=o,w
(3.23)

После подстановки выражения фазой скорости через суммарную скорость

uw = fwuT + kfwλo(∇pc + (ρo - ρw)ge3)
(3.24)

в уравнение сохранения водяной фазы

  ∂sw              uw       sw∂ (ϕρw )
ϕ ∂t-+ ∇uw  = qw - ρw ∇ ρw - ρw--∂t- ,
получим уравнение для вычисления насыщенности:
ϕ ∂sw-+ uT∇fw + ∇ (kfw λo∇pc) = qw - uw-∇ρw - sw-∂(ϕρw)
  ∂t                               ρw       ρw   ∂t
- fw∇uT - ∇ (kfwλo(ρo - ρw )ge3),
или
  ∂sw-     ′               ′
ϕ ∂t  + uTfw∇sw - ∇ (kfwλo|pc|∇sw)
= qw- uw∇ ρw- sw-∂(ϕ-ρw)- fw∇uT - k(ρo- ρw)∇ (fw λo)ge3- fwλo∇(k(ρo- ρw)ge3),
      ρw      ρw   ∂t
или
ϕ ∂sw-+ (uTf ′w + k(ρo - ρw )(fwλo)′ge3)∇sw - ∇ (kfw λo|p′c|∇sw )
   ∂t
= q  - uw-∇ρ  - sw-∂(ϕ-ρw)- f ∇u   - f λ ∇(k(ρ - ρ )ge ).
   w   ρw   w   ρw   ∂t     w   T    w o     o   w   3
(3.25)

Таким образом получили систему уравнений:

  ∂sw
ϕ ----+ (uTf ′w + k(ρo - ρw )(fwλo)′ge3)∇sw - ∇ (kfw λo|p′c|∇sw )
   ∂t
       u        s  ∂(ϕ ρ )
= qw - -w-∇ρw - -w-----w-- fw∇uT  - fwλo∇(k(ρo - ρw)ge3),
       ρw       ρw   ∂t
 ∑   sα∂(ϕρα) ∂p-        ∑    uα-      ∑
α=o,w ρα  ∂p   ∂t + ∇uT + α=o,w ρα∇ ρα = α=o,wqα.
(3.26)

Выберем в качестве расчетных переменных водонасыщенность

s = sw
и глобальное давление
        ∫       ′
p = po -      fw (pc(ξ))dξ  (λT∇p = λw ∇pw + λo∇po).
         pc(s)
Запишем фазовые скорости фильтрации:
uw = - kλw(∇p - fo∇pc + ρwge3),
uo = - kλo(∇p + fw∇pc + ρoge3).
Суммарную скорость можно записать как:
u  = - k(λ ∇p + (λ ρ + λ ρ )ge ).
 T       T       w w    o o  3
(3.27)

Подставим выражение для суммарной скорости в уравнение для давления и получим:

(              )
  ∑   sα-∂(ϕρα)  ∂p- ∇(kλ ∇p )- ∇ (k(λ ρ +λ ρ )ge )+  ∑   uα∇ ρ  = ∑   q
 α=o,wρα   ∂p    ∂t      T          w  w  o o   3  α=o,w ρα   α   α=o,w α
или
( ∑            )                                    ∑
      sα-∂(ϕρα)  ∂p- ∇(kλT∇p )- k(λw∇ ρw+ λo∇ρo)ge3+    uα-∇ρα =
 α=o,wρα   ∂p    ∂t                                α=o,w ρα
 ∑   q + (ρ ∇ (kλ )+ ρ ∇(kλ ))ge .
α=o,w α    w     w    o    o    3
(3.28)

Преобразуем уравнении для насыщенности:

ϕ ∂sw+ ∇ (kf λ ∇p )- ∇ (kλ (∇p + ρ ge )) =
  ∂t        w o  c        w       w  3
       uw       sw ∂(ϕρw)
- qw - ρw∇ ρw - ρw--∂t--.
или
 ∂sw-
ϕ ∂t + ∇ (kfw λo∇pc)- k(∇p + ρwge3)∇ λw =
- qw - uw-∇ρw - sw-∂(ϕ-ρw)+ λw∇ (k∇p) +λw ∇(kρw)ge3.
      ρw       ρw   ∂t

Таким образом система уравнений выглядит следующим образом:

  ∂sw-          ′                     ′
ϕ ∂t  - ∇ (kfwλo|pc|∇sw) - k(∇p + ρwge3)λw∇sw =
- qw - uw-∇ρw - sw-∂(ϕ-ρw)+ λw∇ (k∇p) +λw ∇(kρw)ge3,
      ρw       ρw   ∂t
(              )
  ∑   sα-∂(ϕρα)  ∂p-                                ∑   uα-
      ρα   ∂p    ∂t- ∇(kλT∇p )- k(λw∇ ρw+ λo∇ρo)ge3+    ρα ∇ρα =
 α=o,w                                             α=o,w
   ∑
-      qα + (ρw∇ (kλw)+ ρo∇(kλo))ge3.
  α=o,w
(3.29)

Запишем полученную систему в виде уравнений конвекции-диффузии для насыщенности (s = sw):

as(x,s,p)∂s+ ∇ (ds(x,s,p)∇s)+ cs(x,s,p)∇s = fs(x,s,p),
        ∂t
(3.30)

где

as(x,s,p) = ϕ(x), (3.31)
ds(x,s,p) = - k(x)f w(s)λo(s)|pc(s)|, (3.32)
cs(x,s,p) = - k(x)(p + ρ w(p)ge3)λw(s), (3.33)
fs(x,s,p) = - q w -u
-w-
ρwρw -s
---
ρw∂(ϕ(x)ρ )
-------w-
    ∂t (3.34)
+ λw(s)(k(x)p) + λw(s)(k(x)ρw)ge3. (3.35)
и для давления (p)
         ∂p
ap(x,s,p)--+ ∇ (dp(x,s,p)∇p)+ cp(x,s,p)∇p = fp(x,s,p),
         ∂t
(3.36)

где

ap(x,s,p) = α=o,ws
-α-
ρα∂(ϕρ )
----α-
  ∂p = (scw + soco)ϕ + ϕcr, (3.37)
dp(x,s,p) = - k(x)λ T(s), (3.38)
cp(x,s,p) = (u wcw + uoco) (3.39)
- k(x)(λw(s)ρw(p)cw + λo(s)ρo(p)co)ge3, (3.40)
fp(x,s,p) = - α=o,wqα + (ρw(w) + ρo(o))ge3. (3.41)

Оба уравнения, как для насыщенности так и для давления приводят нас к нестационарным уравнениям конвекции-диффузии в недивергентном виде, следующего вида

         3             3     (         )
a(x)∂u+ ∑  c (x)-∂u-- ∑  -∂-- d (x)-∂u-  = f(x),  x ∈ Ω,  t > 0,
    ∂t  α=1 α   ∂xα   α=1∂xα   α   ∂xα
(3.42)

при стандартных предположениях 0 < k1 k(x) k2.

С граничными условиями Неймана

    ∂u(x,t)-
d(x )  ∂n   = 0, x ∈ δΩ,  t > 0.
(3.43)

Для однозначной разрешимости нестационарной задачи задается начальное условие

u(x,0) = u0(x), x ∈ Ω.
(3.44)

На множестве функций u(x), удовлетворяющих граничным условиям (3.43), нестационарную задачу конвекции-диффузии запишем в виде дифференциально-операторного уравнения

∂u+ Au = f,  A = C + D,  t > 0.
∂t
(3.45)

Здесь D - оператор диффузионного переноса, который определяется выражением

        ∑3 -∂--(     -∂u-)
Du  = -    ∂xα  dα(x)∂xα  .
       α=1
Для оператора конвективного переноса в недивергентной форме имеем
          ∑3
C = C1u =    cα(x ) ∂u-.
          α=1     ∂xα

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

Обозначим через ω равномерную сетку в области Ω

ω = {x|x = (x1,x2,x3), xα = iαh α, iα = 0,1,...,Nα,  N αhα = lα}
причем ω — множество внутренних узлов, а δω - множество граничных узлов.

Запишем разностный оператор диффузионного переноса

     3                 (       )
D = ∑   D  ,  D   y = -  a(α)yx    ,  α = 1,2,3,  x ∈ ω,
    α=1  (α)    (α)            α xα
где для гладких коэффициентов диффузии можно положить
a(1)(x) = k(x1 - 0.5h1,x2,x3),
a(2)(x) = k(x1,x2 - 0.5h2,x3),
a(3)(x) = k(x1,x2,x3 - 0.5h3).
Заметим, что данный дискретный аналог аппроксимирует дифференциальный оператор с точностью O(h2).

Для недивергентного оператора конвективного переноса можно использовать аппроксимации центральными разностными производными

      ∑3
C1y =    C(1α), C (α1)y = b(α)y˙xα, α = 1,2,3, b(α)(x) = vα(x ),  x ∈ ω.
      α=1
Такой разностный оператор конвективного переноса имеет второй порядок аппроксимации.

Более лучшими являются аппроксимации оператора конвективного переноса, которые основаны на использовании аппроксимаций с заданием коэффициентов vα(x) на сдвинутых на полшага сетках в соответствующих направлениях.

Будем относить коэффициент v1(x) конвективного переноса по переменной x1 к узлам сетки, которая сдвинута по этой переменной на полшага. Сетка для коэффициента v2(x) сдвигается по x2 на 0.5h2. Аналогично сетка сдвигается и для третьего направления. Интегро-интерполяционным методом для разностного оператора конвективного переноса в недивергентной форме получим

 (1)  1 ( (1)                     (1)                   )
C1  = 2  b  (x1 - 0.5h1,x2,x3)yx1 + b (x1 + 0.5h1,x2,x3)yx1 ,
      1 (                                             )
C(22)= -  b(2)(x1,x2 - 0.5h2,x3)yx2 + b(2)(x1,x2 + 0.5h2,x3)yx2 ,
      2
        (                                             )
C(2)= 1  b(3)(x1,x2,x3 - 0.5h3)yx3 + b(3)(x1,x2,x3 + 0.5h3)yx3 ,
 3    2
     ∑3  (α)
C1 =    C1  , x ∈ ω.
     α=1

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

bα(x) = b+(x)+ b-(x),
        α      α
 +      1
bα(x) = 2 (bα(x)+ |bα(x)|) ≥ 0,
       1
b-α(x) =- (bα(x)- |bα(x)|)) ≤ 0.
       2
поставим в соответсвие дифференциональному оператору простейший разностный оператор конвективного переноса C = C1 с направленными разностями
      ∑3
C1y =    (b+α(x)yxα + b-α(x)yxα), x ∈ ω.
      α=1
который имеет только первый порядок аппроксимации.

Таким образом от уравнения (3.45) придем к дифференциально-операторному уравнению

dy+ Ay = ϕ,  A = C + D,  0 < t < T,
dt
(3.46)

на множестве сеточных функций y(t) H с начальным условием

y(t) = y0.
(3.47)

Для численного решения уравнения 3.46 применяются схемы с весами вида:

yn+1 --yn+C (σ1yn+1+ (1- σ1)yn)+D (σ2yn+1+(1- σ2)yn) = ϕn, tn ∈ ω τ
    τ
y  = u .
 0    0
(3.48)

3.3.3 Реализация

Для задания параметров присущих данной математической модели создадим класс Parametrs

// Parameters.h 
class Parameters { 
public: 
    Parameters(); 
 
    double se(double sw); 
    double pcapillary(double sw); 
    double dpc(double sw); 
    double krw(double sw); 
    double krn(double sw); 
 
    double muWater, muOil; 
    double psi, cWater, cOil, cRock; 
    double s_rWater, s_rOil, rhoWater, rhoOil; 
    double g; 
    double teta, pe; 
}; 
    

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

// Parameters.cpp 
Parameters::Parameters(){ 
    muWater = 0.001; 
    muOil   = 0.01; 
 
    psi     = 6894.757; 
    cWater  = (3 * 1.e-6 / psi ); 
    cOil    = (22.14 * 1.e-6 / psi ); 
    cRock   = (6 * 1.e-6 / psi ); 
 
    s_rWater  = 0.2; 
    s_rOil    = 0.15; 
    rhoWater  = 1000.0; 
    rhoOil    = 900.0; 
    g         = 9.8; 
 
    teta = 2; 
    pe   = 1.0e4; 
} 
    

Для определения коэффициентов зависимых от насыщенности применяется модель Брукс-Корей:

double Parameters::se(double sw) { 
    if (sw <= s_rWater) 
        return 0.; 
    if (sw >= (1 - s_rOil)) 
        return 1.; 
    return (sw - s_rWater) / (1 - s_rWater - s_rOil); 
} 
 
double Parameters::pcapillary(double sw) { 
    return pe * pow(se(sw), -1.0 / teta); 
} 
 
double Parameters::dpc(double sw) { 
    if (sw <= s_rWater) 
        sw = s_rWater + 1.e-5; 
    if (sw >= (1.0 - s_rOil)) 
        sw = 1.0 - s_rOil - 1.e-5; 
 
    double v = 1.0 / teta * pe * pow(se(sw), (-1.0 / teta - 1)) * 1.0 / (1.0 - s_rOil - s_rWater); 
    return v; 
} 
 
double Parameters::krw(double sw) { 
    double lw; 
    lw = pow(se(sw), (2 + 3 * teta) / teta); 
    return lw; 
} 
 
double Parameters::krn(double sw) { 
    double ln; 
    ln = pow(1 - se(sw), 2) * (1 - pow(se(sw), (2 + teta) / teta)); 
    return ln; 
} 
    

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

// PressureDiffOperator.h 
class PressureDiffOperator : public DiffusionOperator { 
public: 
    PressureDiffOperator(string prevPressure, Variable& saturation, 
                                Variable& permX, Variable& permY, Variable& permZ, 
                                Variable& wellWater, Variable& wellOil); 
    virtual ~PressureDiffOperator(); 
 
    virtual void prepare(int level); 
 
protected: 
    virtual double diff(int i, int j, int k, int di, int dj, int dk); 
 
private: 
    string mSaturationName; 
    string mPressureName; 
    DAValue* mSaturation; 
    DAValue* mPressure; 
 
    string mPermXName; 
    string mPermYName; 
    string mPermZName; 
    string mWellOilName; 
    string mWellWaterName; 
    DAValue* mPermX; 
    DAValue* mPermY; 
    DAValue* mPermZ; 
    DAValue* mWellOil; 
    DAValue* mWellWater; 
 
}; 
    

Реализации конструктора класса PressureDiffOperator принимает все зависимые переменные

// PressureDiffOperator.cpp 
PressureDiffOperator::PressureDiffOperator(string prevPressure, Variable& saturation, 
                                Variable& permX, Variable& permY, Variable& permZ, 
                                Variable& wellWater, Variable& wellOil) 
        : DiffusionOperator() 
        , mSaturationName(saturation.getName()) 
        , mPressureName(prevPressure) 
        , mPermXName(permX.getName()) 
        , mPermYName(permY.getName()) 
        , mPermZName(permZ.getName()) 
        , mWellWaterName(wellWater.getName()) 
        , mWellOilName(wellOil.getName()){ 
} 
 
PressureDiffOperator::~PressureDiffOperator() { 
} 
    

Метод prepare в классе PressureDiffOperator вытаскивает ссылки из Storage для зависимых переменных

// PressureDiffOperator.cpp 
void PressureDiffOperator::prepare(int level){ 
    DiffusionOperator::prepare(level); 
 
    mSaturation = &(Storage::instance().get(mSaturationName).getDAValue(level)); 
    mPressure   = &(Storage::instance().get(mPressureName).getDAValue(level)); 
 
    mPermX   = &(Storage::instance().get(mPermXName).getDAValue(level)); 
    mPermY   = &(Storage::instance().get(mPermYName).getDAValue(level)); 
    mPermZ   = &(Storage::instance().get(mPermZName).getDAValue(level)); 
 
    mWellOil = &(Storage::instance().get(mWellOilName).getDAValue(level)); 
    mWellWater = &(Storage::instance().get(mWellWaterName).getDAValue(level)); 
} 
    

В реализации класса PressureDiffOperator задаем коэффициент диффузии с учетом граничных условий Неймана

// PressureDiffOperator.cpp 
double PressureDiffOperator::diff(int i, int j, int k, int di, int dj, int dk){ 
    if(   (i == 0        && di == -1) 
       || (i == (mN.x-1) && di ==  1) 
       || (j == 0        && dj == -1) 
       || (j == (mN.y-1) && dj ==  1) 
       || (k == 0        && dk == -1) 
       || (k == (mN.z-1) && dk ==  1) 
       ) 
        return 0.0; 
    //Newman boundary condition 
 
    DAValue* permXYZ; 
    double h; 
    if(di != 0){ 
        permXYZ = mPermX; 
        h = mH.x; 
    } 
    if(dj != 0){ 
        permXYZ = mPermY; 
        h = mH.y; 
    } 
    if(dk != 0){ 
        permXYZ = mPermZ; 
        h = mH.z; 
    } 
 
    Parameters p; 
    double perm = (permXYZ->get(i + di, j + dj, k + dk)  + permXYZ->get(i, j, k)) / 2; 
    double s = (mSaturation->get(i + di, j + dj, k + dk) + mSaturation->get(i, j, k)) / 2; 
    double lambdaWater = p.krw(s) / p.muWater; 
    double lambdaOil   = p.krn(s) / p.muOil; 
 
    return -(lambdaWater + lambdaOil) * perm; 
} 
    

Разностный оператор конвекции для давления реализуется посредством наследования класса ConvectionOperator:

// PressureConvOperator.h 
class PressureConvOperator : public ConvectionOperator { 
public: 
    /* set pressure=oldPressure becouse we solve usong linear solver*/ 
    PressureConvOperator(Approximation app, string prevPressure, Variable& saturation, 
                                Variable& permX, Variable& permY, Variable& permZ, 
                                Variable& wellWater, Variable& wellOil); 
    virtual ~PressureConvOperator(); 
 
    virtual void prepare(int level); 
 
protected: 
    virtual double conv(int i, int j, int k, int di, int dj, int dk); 
 
    private: 
    string mSaturationName; 
    string mPressureName; 
    DAValue* mSaturation; 
    DAValue* mPressure; 
 
    string mPermXName; 
    string mPermYName; 
    string mPermZName; 
    string mWellOilName; 
    string mWellWaterName; 
    DAValue* mPermX; 
    DAValue* mPermY; 
    DAValue* mPermZ; 
    DAValue* mWellOil; 
    DAValue* mWellWater; 
}; 
    

Реализация оператора конвективного переноса для уравнения давления аналогична оператору диффузии