Численные алгоритмы классической матфизики. XI. О вычислении собственных значений уравнения переноса

16 Pages·0.233 MB·Russian

by  Алгазин С.Д.

Институт проблем механики
Российской Академии Наук

С. Д. Алгазин

ЧИСЛЕННЫЕ АЛГОРИТМЫ КЛАССИЧЕСКОЙ
МАТФИЗИКИ.

XI. О вычислении собственных значений уравнения
переноса.

Препринт № 801

Москва 2006 г.
Аннотация.

Рассматривается задача на собственные значения для
уравнения второго порядка с переменными коэффициентами.
Построен численный алгоритм без насыщения. Приводится текст
программы на Фортране-77.

The abstract.

The problem on eigenvalues for a second-kind equation with
floating factors is considered. The numerical algorithm without saturation
builts. The text of the program on the FORTRAN - 77 is reduced.

Работа выполнена при финансовой поддержке Российского
фонда фундаментальных исследований. Проект № 05-01-00250.

055(02)2  Институт проблем механики РАН 2006

2
Введение.

В [1] рассматриваются задачи на собственные значения для оператора
Лапласа в произвольной гладкой области с постоянными коэффициентами.
Однако ряд задач математической физики приводит к задачам на собственные
значения для уравнения второго порядка с переменными коэффициентами (см.
ниже). Для решения этих задач существует метод наискорейшего спуска [2],
который в частности сводит решение самосопряжённого уравнения второго
порядка к последовательности решения задач для уравнения Пуассона в этой же
области. Этот метод применяется также для решения нелинейных уравнений [3].
Однако рассмотренные там примеры численных расчётов не вызывают
оптимизма в быстроте сходимости метода. В настоящей главе построен
численный алгоритм без насыщения в задаче на собственные значения для
эллиптического уравнения второго порядка с переменными коэффициентами.
Для примера рассмотрено краевое условие Неймана. По ходу изложения будет
пояснено, как рассмотреть другие краевые условия.

§ 1. Постановка задачи фильтрации газа в пористой среде.

Искомое уравнение имеет вид:

∂(mρ) 
(1.1) +div(ρv) =0,
∂t

где m=V /V – пористость ( для реальных пластов лежит в пределах 0,15 ~ 0,22
пор.

); mρ – концентрация; v - скорость фильтрации (а не скорость жидкости).
Это уравнение получается из обычного закона сохранения массы

d d
(1.2) ∫ρdτ= ∫ρmdτ=0,
dt dt
V V
пор.

где V – объём пор, а V – полный объём, причём оба объёма подвижные. Из
пор.
(1.2) получаем, применяя формулу дифференцирования по подвижному объёму
[4]:

∂(mρ)   
= div(mρw), v = mw,
∂t

 
где v − скорость фильтрации, а w−скорость жидкости. В результате получаем
уравнение (1.1).
Закон Дарси (1856) справедлив для медленных движений жидкости в
изотропной пористой среде, т.е. для малых значений числа Рейнольдса Re
(Re<Re )
кр.

 k
(1.3) v =− qrad p,
µ
3
где k- коэффициент проницаемости, измеряемый в Дарси (1д = 10-8/0,981 см2).
Для реальных пористых сред k = 100 ~ 1000 мд ( 1мд=10-3д). Проницаемость -
геометрическая характеристика пористой среды, т. е. определяется размерами
частиц, их формой и упаковкой, μ – динамическая вязкость.
Уравнение состояния.

M p
ρ= ,
RT z(p)

где M – молярный вес газа, R – универсальная газовая постоянная, T –
абсолютная температура; z(p) – определяется экспериментально (z(p) =1 для
совершенного газа), т.е. это баротропный газ.
Уравнение (1.1) относится к случаю, когда в пласте нет источников газа
(скважин). В общем случае уравнение неразрывности имеет вид:

∂(mρ) 
(1.4) +div(ρv)= f(z,t),z∈G,
∂t

где f(z,t) – заданная функция, G – двумерная область с гладкой границей дG∈C∞.
Пусть z =ϕ(ς),ς= reiθ - конформное отображение единичного круга на область
G. Выпишем уравнение (1.5) в новых переменных [4]:

ds2 =(dr2 +r2dθ2)|ϕ′(ς)|2⇒ g =|ϕ′(ς)|2, g = r2 |ϕ′(ς)|2, g =|ϕ′(ς)|2 r.
11 22
1 ∂p
grad p = ,
r |ϕ′(ς)| ∂r

1 ∂p
gradp = .
θ |ϕ′(ς)|r ∂θ

Подставляя в (1.4) получим

∂(mρ)
(1.5) =|ϕ′(ς)|−2 L(w)+ f(ς,t),ς= reiθ,0≤ r ≤1,0≤θ< 2π,|ς|≤1;
∂t

1 ∂  ∂w 1 ∂  ∂w
(1.6) L(w) = rk(r,θ) + k(r,θ) .
r ∂r ∂r  r2 ∂θ ∂θ

Здесь введены обозначения:

m=m(r,θ) – пористость (известная функция координат );
p=p(r,θ,t) – давление ( неизвестная функция координат и времени );
k =k(r,θ,p)=k(r,θ)ψ(p)- проницаемость ( известная функция координат и
давления);
ρ=ρ(p) – плотность ( известная функция давления );
4
μ=μ(p) – вязкость ( известная функция давления );
ρ(p)ψ(p)
w(p) = ∫ dp.
µ(p)

Размерность: (M – единица массы, L – единица длины, T – единица
времени).

m,θ,ψ – безразмерные величины;
[p]=M/LT2, [ρ]=M/L3, [μ]=M/LT, [k]=L2, [w]=M/L3T, [r]=L
f(ς,t)=f(r,θ,t) – плотность отбора газа, т.е. масса газа, выделяющаяся в единицу
времени в единице объёма в пласте. Если ввести мощность пласта h=h(x,y) (т. е.
высоту пласта в точке (x,y)∈G, рассматриваемой области), то вид уравнения
(1.5) не изменится, если заменить m на mh, а k на kh. В этом последнем случае
[f]=M/L2T, т.е. масса, выделяющаяся из пласта в единицу времени и с единицы
площади.
Таким образом, (1.5), (1.6) искомая постановка задачи фильтрации. К
этому уравнению нужно добавить граничное условие

∂p
(1.7) =0,
∂n
∂G
которое означает отсутствие потока газа через границу области дG (см. (1.3)).
Заметим, что функция w также удовлетворяет этому граничному условию.

§2. Дискретизация по пространственным переменным.

Для дискретизации задачи (1.5) – (1.7) проведём вначале дискретизацию
оператора L(w). Рассмотрим спектральную задачу:

∂w
(2.1) L(w)+λw=0, =0.
∂r
r=1

Заметим, что

 ∂w2 k ∂w2
− ∫ L(w)wdς= ∫ k  +   dς.
|ς|≤1 |ς|≤1  ∂r  r2 ∂θ 

Таким образом, краевая задача (2.1) эквивалентна следующей
экстремальной задаче

 ∂w2 k ∂w2 
(2.2) J(w) = ∫ k  +   −λw2dς →min.
|ς|≤1  ∂r  r2 ∂θ 

Действительно, δJ (вариация функционала J) есть главная линейная часть
приращения J(w+h)-J(w), где h – произвольная гладкая функция. Нетрудно
получить, что
5
 k  1 ∂ 1 ∂
δJ = 2 ∫ kw h + w h −λwh dς= 2{krw h − ∫ [ (rkw )+ (kw )+
 r r r2 θ θ  r r=1 r ∂r r r2 ∂θ θ
|ς|≤1 |ς|≤1
+λw]hdς}=0

Так как h – произвольная функция, отсюда следуют соотношения (2.1). Итак,
при поиске минимума функционала (2.2) не нужно заранее удовлетворять
краевому условию Неймана, т. е. это краевое условие естественное.
Для дискретизации функционала (2.2) применим квадратурную формулу
[1]:

∫ f(ς)dσ= ∑c f , f = f(r eiθl),
νl νl νl ν
|ς|≤1 ν,l
(2.3)
(2ν−1)π 2πl
r =cos ,ν=1,2,...,m;θ = , l =0,1,...,2n; N = 2n+1.
ν 4m l N
Эта квадратурная формула получается, если заменить подынтегральную
функцию интерполяционной формулой для функции двух переменных в круге
[1]:
2n m
(2.4) (P f)(r,θ) = ∑∑ f L (r,θ), f = f(r ,θ),
M νl νl νl ν l
l=0ν=1
2T (r) D (θ−θ) D (θ−θ +π)
L (r,θ) = 2m  n l − n l ,
vl NT' (r ) r −r r +r 
2m ν ν ν
n
D (θ)=0,5+∑coskθ, T (r) =cos(marccosx).
n m
k=1
Интерполяционная формула (2.4) обладает нужными свойствами.
Действительно, формула (2.4) точна на многочленах от двух переменных
степени ω=min(n,m-1). Обозначим множество этих многочленов P , а E
ω ω
обозначим наилучшее приближение функции f∈C[D] (D – единичный круг)
многочленом из P . Тогда определён проектор
ω
P : C[D]→LM, LM=L(L ,…, L )
M 1 M
и справедливо классическое неравенство:
(2.5) | f(r,θ)−(P f)(r,θ)|≤(1+| P | )E (f),
M M ∞ ω
в котором |P | - норма проектора P . Так же, как в одномерном случае,
M∞ M
неравенство (2.5) показывает, что соответствующая интерполяционная формула
не имеет насыщения. Норма проектора P удовлетворяет соотношению
M
|P | =O(ln2M),
M∞
причём не составляет труда уточнить эту оценку. Делая некоторые
предположения о гладкости класса интерполируемых функций, можно оценить
скорость убывания наилучшего приближения E при M→∞ и получить
ω
конкретные оценки погрешности интерполяционной формулы (2.4).
Пусть
f(r,θ)=(P f)(r,θ) +ρ (r,θ;f),
M M
6
где ρ (r,θ;f) - погрешность интерполяционной формулы (2.4) (остаток). Тогда
M
справедлива следующая теорема К. И. Бабенко.
Теорема (К. И. Бабенко). Рассмотрим класс функций
HM(K;D)⊂C(D),удовлетворяющих в круге D условиям
∞
∂k+l f
≤ K, k +l ≤µ,
∂xk∂yl
тогда, если f∈HM(K;D), то
∞
(2.6) |ρ ( . ;f)| ≤ c K M−µ/2 log2M,
M ∞ µ
где c - константа, зависящая от µ.
µ
Таким образом, из рассмотрения формулы (2.6) видно, что при
одинаковом числе узлов интерполяции M скорость убывания погрешности
интерполяционной формулы (2.4) возрастает с ростом µ, т.е. с ростом гладкости
интерполируемой функции f. Это означает, что полученная интерполяционная
формула не имеет насыщения.
Основываясь на интерполяционной формуле (2.4), легко построить
квадратурную формулу для вычисления определённых интегралов, когда
областью интегрирования является круг. В самом деле, заменяя
подынтегральную функцию выражением (2.4), получим квадратурную формулу
(2.3), где dσ - элемент площади, с – весовые коэффициенты, а δ(f) –
νl
погрешность. Для с имеем
νl
c = ∫L (r,θ)dσ,
νl νl
D
и они не зависят от l. Введём в рассмотрение блочно-диагональную матрицу
C=diag(c , c ,…, c ),
1 2 m
где c , ν=1,2,…,m – диагональные матрицы размера N×N с одинаковыми
ν
числами на диагонали. Для погрешности квадратурной формулы имеем
следующую оценку:
|δ(f)|≤ 2πE (f).
ω
Заметим, что все c положительны при достаточно большом числе узлов
νl
интерполяции.
Для коэффициентов квадратурной формулы (2.3) имеем выражение:
4πr cosψ m−1  s−1 (2ν−1)π
с = ν  ν + ∑t cossψ , t =1/(1+(−1) 2 s),ψ = , s≥1 −нечётно.
ν m(2nν +1) 2 s=3(2)s ν s ν 4m

∂w
  = ∑H w
 ∂r  νl,µp µp
ς=ς µ,p
νl

∂w N
  = ∑B w
∂θ lp νp
ς=ς p=1
νl

7
Матрицы B и H подучаются дифференцированием интерполяционной
формулы (2.4)

2 n 2π(l − p)
B = ∑ksink
lp N N
k=1

Для получения матрицы H продифференцируем по r интерполяционную
формулу (2.4). Обозначим
d  T (r)  1 2m−1scossψ sinsψ
A(1) =  2m  = ∑ ν µ
µν dr(r −rν)T2′m(rν)r=rµ m s=1 sinψµ
d  T (r)  1 2m−1s(−1)s cossψ sinsψ
A(2) =  2m  = − ∑ ν µ
µν dr(r +rν)T2′m(rν)r=rµ m s=1 sinψµ
Дифференцируя (2.4) по r получим
du(r,θ) m  2 2n 
= ∑A(1)u − ∑A(2)D (θ +π−θ)u 
dr θr==θrµ ν=1 µν νp N l=0 µν n p l νl
p
n
где D (θ +π−θ) =0,5+∑(−1)k cosk(θ −θ) ⇒
n p l p l
k=1
2
H = A(1)δ − A(2)D (θ +π−θ)
µp,νl µν pl N µν n p l
Нетрудно видеть, что H – h-матрица [1], и, следовательно, представляется в
виде
2 n '
(2.7) H = ∑ Λ ⊗h ,
N k k
k=0
где штрих у знака суммы означает, что слагаемое при k=0 берётся с
коэффициентом ½, Λ, k=0,1,…,n – матрица размера m × m; h , k=0,1,…,n –
k k
матрица размера N × N:

h =cos[k2π(i-j)/N)], i,j=1,2,…,N,
kij

через ⊗ обозначено кронекерово произведение матриц. Конкретно матрицы Λ в
k
этом случае имеют вид:

Λ =(−1)k+1A(2) + A(1)
kµν µν µν
Итак,

2 2m−1 scossψ sinsψ
Λ = ∑ ν µ
2k,µν m sinψ
s=2(2) µ
8
2 2m−1 scossψ sinsψ
Λ = ∑ ν µ
2k+1,µν m sinψ
s=1(2) µ
Ниже будем обозначать эти матрицы - Λ(k). Распишем формулу (2.7)
µν
подробно:
2 n ' 2π(p−l)
H = ∑ Λ(k) cosk ,
νl,µp N νµ N
k=0
~
2 n ' 2π(l −l )
H = ∑ Λ(q) cosq .
νl,µ~~l N νµ~ N
q=0
Используя квадратурную формулу (2.3) функционал (2.2) преобразуем в
квадратичную форму:
 ∂w2 k ∂w2 
(2.8) J(w) = ∑c k   + νl   +λw2.
ν,l νl νl ∂r ς=ς rν2 ∂θς=ς νl
νl νl
Дифференцируя (2.8) по w получим
µ~~l
N
∑B w +∑Aµ~w −λc w =0,
µ~~l,µp µp ~lp µ~p µ~ µ~~l
µ,p p=1
где
~
4 n ' n ' m 2n 2π(p−l) 2π(l −l )
Bµ~~l,µp = N2 ∑k=0 ∑q=0 ν∑=1cνΛ(νkµ)Λ(νqµ~)∑l=0kνl cosk N cosq N 
~
A~lµ~p = N42 crµµ~~2 ∑kn=1∑qn=1kq∑lN=1 kµ~l sink 2π(lN− p)sinq2π(Nl −l ).
Это есть дискретный аналог задачи на собственные значения
div(kqradw)+λw=0, r<1
∂w
=0.
∂r
r=1
Оценка погрешности описанной дискретизации может быть получена по
схеме, описанной в [1]. См. также [5].
§3. Результаты численных экспериментов.
Проводились расчёты задачи (2.1) в круге (k=1), возмущённом круге
(эпитрохоида, k=1, k≠1). В круге собственные значения известны λ,i =1,2,...-
i
нуль производной функции Бесселя. Сравнение с таблицами показывает, что
даже на сетке 3×7 имеем в первых собственных значениях по 4 знака после
запятой. Однако точность значительно хуже, чем по методике, описанной в [1]
для уравнения с постоянными коэффициентами. Второй расчёт проводился для
9
эпитрохоиды (ϕ(ς)=ς(1+εςn), ε = 0.0625, n=12), для которой в [1] в таблице (3.9)
приведены вычисленные значения. Результаты расчётов представлены в
таблице 3.1 (выписаны знаки, совпавшие с расчётами из [1]).

Таблица 3.1
№ 8×11 10×21 15×31
2 1.76 1.7751 1.77623557
6 3.77 3.72 3.7317
11 5.2 5.16 5.1770
16 6.9 6.5 6.4787
Третий расчёт проводился для той же эпитрохоиды для функции
k(r,θ)=k (0.1+r2)(sin12θ+1.1), k =10-13/0.981 м2=0.1 дарси.
0 0
Результаты расчётов представлены в таблице 3.2.
Таблица 3.2
№ 8×11 10×21 15×31 20×41 30×41
2 2.8451E-7 2.5447E-7 2.4840E-7 2.6227E-7 2.5541E-7
6 5.8268E-7 6.5814E-7 7.0164E-7 7.2177E-7 7.1376E-7
11 9.2409E-7 9.3812E-7 1.0838E-6 9.8923E-7 9.5191E-7
16 1.2075E-6 1.3895E-6 1.1905E-6 1.1763E-6 1.0876E-6
21 1.4927E-6 1.5927E-6 1.4384E-6 1.4482E-6 1.4222E-6

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

§4. Описание программного комплекса.

Решение поставленной задачи осуществляет программа ALAPN:

$objcomment lib:"treq.lib"
PROGRAM ALAPN
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION H(1512900),B(1512900),WR(1230),WI(1230),IANA(1230),
*U(1230),PR(1230),A1(900),A2(900),NL(30),Z(1230),X(1230),
*Y(1230),Y2(11)
REAL*4 PR
P(R,T)=(1.D-13/0.981)*(0.1+R*R)*(SIN(12.*T)+1.1)
C P(R,T)=1.D0
PI=3.141592653589D0
WRITE (*,*) 'M=?, N=?'
READ (*,*) M,N
10