Алгоритм приближенного решения задачи оптимального управления уровнем грунтовых вод при слоистом строении водоносных пластов

Бюллетень науки и практики / Bulletin of Science and Practice

Рассмотрим движение подземных вод в двухслойных водоносных пластах, состоящих из покровной толщи, подстилаемой основным напорным водоносным горизонтом, разделенным от нижележащего напорного пласта слабопроницаемой прослойкой. Уровни грунтовых вод (УГВ) в покровном слое и напорных вод в первом напорном пласте обозначим соответственно h(x, y, t) и H(x, y, t), а коэффициенты фильтрации и мощности этих пластов - соответственно kb(x,y), mb(x,y, t) = h(x, y, t) — b(x,y) и k(x,y), m(x,y), где b(x, y) - поверхность раздела между покровным слоем и первым сверху напорным слоем.

Течение подземных вод в таких слоистых пластах описывается следующей системой дифференциальных уравнений [1].

_μb ∂h_- ∂_(Tb∂h_)-∂_(Tb∂h_)+kb h-H _=fb, ∂t ∂x ∂x ∂y ∂y       m_b

∂H ∂   ∂H ∂   ∂H

μупр ∂t - ∂x ( ∂x) - ∂y ( ∂y

-

k b

h-H m b

k

+ ⁿ (H - Z) = f, m n

(x, y) ∈ D,t > 0

с начально-краевыми условиями

h(x, y, 0) = h 0 (x, y),

H(x, y, 0) = H 0 (x, y), (x, y) ∈ d,

∂h

Tb   + βbh = αb,

T + βHh = α, (x, y) ∈ S = ∂D,t > 0, ∂n где μb и μупр — свободная водоотдача в покровном слое и упругая водоотдача в первом напорном горизонте;

T b (x, y, t) = k b (x, y) ∙ m b (x, y, t) = k b (h -b)                                 (4)

T(x, y) = k(x, y) ∙ m(x, y) - водопроводимость этих слоев; h₀(x,y) и H₀(x,y) — начальное распределение УГВ и напоров; Z(x, y) — напоры во втором напорном пласте; k_m , m_n — коэффициент фильтрации и мощность слабопроницаемой прослойки между напорными пластами; f_b(x, y, t) — функция инфильтрации; f(x, y, t) — функция, учитывающая работу скважин, пробуренных в первый напорный пласт; α_b(x, y, t), β_b(x, y, t), α(x, y, t) и β(x, y, t) — известные функции; D — плоская область в плане, S = ∂D — ее граница; ⃗ n ⃗ — вектор внешней нормали к границе области. Задача оптимального управления УГВ заключается в нахождении функции f_b(x, y, t) , доставляющей при t ≥ T₀ минимум функционалу

J(f b ) = ∬ _D[h(x, y, T 0 , f b (x, y, T 0 )) - φ(x, y)]² dx dy + ∫₀^T0 ∬ _D[f(x, y, t)]² dx dy dt.          ⁽⁵⁾

Здесь h(x, y,t) = h(x, y, t, fb) — УГВ, определяемые из задачи (1)-(3); φ(x, y) — заданная функция, равная оптимальному УГВ; γ > 0 — параметр регуляризации; T0 — заданный момент времени.

Функция h(x, y, t, fb) называется объектом управления, а f(x, y, t) — функцией управления или управлением. Функция fb∗(x, y, t), доставляющая минимум функционалу (5), называется оптимальным управлением, а соответствующая ей функция h∗(x, y, t) = h(x, y, t, fb∗) — оптимальным УГВ. В работе [2] установлено, что оптимальное управление fb∗ должно удовлетворять условию fb(x, y, t) = 1 ψ(x, y, t) ,                                                 (6)

где ψ(x, y, t) - решение сопряженной начально-краевой задачи

μ∂ψ+∂(T∂ψ)+∂(T∂ψ)-∂t ∂x ∂x ∂y ∂y

μ ^∂ζ +^∂(T ^∂ψ )+^∂(T ^∂ψ )+k

⎪ ^упр ∂t ∂x ∂x ∂y ∂y

⎩               (x, y) ∈ D,0 ≤t≤T 0 ,

ψ-ζ kb      = 0, b mb

^ψ-ζ-^kiζ=0, m b    m i

ψ(x, y, T 0 ) = -   [h(x, y, T 0 - φ(x, y)],

μb

⎪ ζ(x, y, T 0 ) = -     [H(x, y, T 0 - H 0 (x, y)], (x, y) ∈ D,

⎩               μ упр

∂ψ

T b   +β b ψ=0,

T ^∂ζ +βζ=0,(x,y) ∈ S,0 ≤t≤T 0 .

Отсюда следует, что при фиксированном fb нужно решить две краевые задачи; сначала из (1)-(3) надо определить функции h(x, y, t, fb) , H(x, y, t) , затем в начальные условия (8) подставить получившиеся функции h(x, y, T0) и H(x, y, T0) и из сопряженной краевой задачи (7)-(9) найти ψ(x, y, t) и полученное ψ(x, y, t) подставить в формулу (6).

Теперь рассмотрим более подробно алгоритмы решения этих задач. Задачи (1)-(3) и (7)(9) решаем методом конечных элементов [3, 4]. Для этого область фильтрации D произвольным образом разбиваем на треугольные элементы и для каждого элемента (е) введем линейные базисные функции

N_s^(e)(x, y) = a s + b s x+ c s y,s = i,j,k

Ns e) (x, y) = as + bsx + csy,     s = i, j, k

, где i,j,k — номера вершин элемента, ai = (xj yk - xkyj)/∆e  , bi = (yi - yk)/∆e  ,  ci = (xk - xj)/∆e  и т.д. Остальные коэффициенты получаются с помощью круговой подстановки индексов i,j,k;

	1   x i    y i
∆ e =	1   x j     y j	- удвоенная площадь треугольника (е)
	1  x k   y k

Искомые функции h(x, y, t) и H(x, y, t) внутри элемента (е) аппроксимируем функциями

h^(e)(x,y,t) = h i (t)N_i^(e)(x, y) + h j (t)N_j^(e)(x, y) + h k (t)N_k^(e) (x, y),                      ⁽¹⁰⁾

H^(e)(x, y, t) = H i (t)N_i^(e)(x, y) + H j (t)N_j^(e)(x, y) + H k (t)N_k^(e)(x, y),

Здесь h s (t) = h(x s ,y s ,t), H s (t) = H(x s , y s , t), s = i, j, k.

Суммируя равенства (10) по всем элементам, получаем аппроксимации искомых функций по всей области D:

h n (x, y, t) = ∑ e ^m =1 h^(e)(x, y, t) = ∑ j ⁿ =1 h j (t)N j (x, y),                            (11)

Hn(x,y,t) =∑em=1H(e)(x,y,t)=∑jn=1Hj(t)Nj(x,y), где m — число всех элементов, n — число всех узлов сетки.

Бюллетень науки и практики / Bulletin of Science and Practice Т. 8. №6. 2022

Далее временной отрезок [0, T0] разобьем на равные промежутки длиной ∆t, так что ts = s∆t, s = 0,1,2,… . На каждом временном слое [ts-1, ts] , s = 0,1,2,… , к системе (1)-(3)

применяем неявную конечно-разностную схему с весом σ (0<σ≤1) [5]:

^μ упр   _s

∆t H

μbhs-σ[∂(T∂hs)+∂(T∂hs)]

∆t       ∂x b ∂x ∂y b ∂y

ks - Hs + σk b mb

= Q^s b ,

∂   ∂Hs ∂   ∂Hs

σ ∂x ∂x    ∂y ∂y

- σk ks - Hs + σ kn (Hs - Z) = Qs, b  mb     mn

σ(T b ^∂hs+β^s b h s ) =A^s b , ∂H^s

σ(T∂n +βsHs) =As,

где h^s = h(x, y, t s ), H^s = H(x,y,t s ), f_b^s = f_b^s(x, y, t s ), f^s = f(x,y,t s ),

Q^s b =σf b ^s+ ^μ ∆ ^b t h s-1 +(1-σ)[ ∂ ^∂ x (T b ^∂h ∂ ^s x ^-1)+ ∂ ^∂ y (T b ^∂h ∂ ^s y ^-1)-k b ^hs-1 m ^- _b ^Hs-1

_{+ fb} s-1

],

Qs=σfs+μупрНs-1+(1-σ)[∂(T∂Hs-1)+ ∂(T∂Hs-1)+k hs-1-Hs-1 ∆t                   ∂x     ∂x      ∂y     ∂y       b    mb ki (Hs-1 mi

-

Z) + fs-1],

A^s_b = σα^s_b + (1 - σ) (α^s_b^-1 - T b ^∂hs-1 - β^s_b^-1h^s-1),

A^s = σα^s + (1 -σ) (α^s -T ^∂hs-1 - β^s-1h^s-1).

В системе (12), (13) вместо h(x, y, t), H(x, y, t) подставим их аппроксимации hn(x, y, t) и Hn(x, y, t) соответственно и применяем обобщенный принцип Галеркина (для удобства записи верхний индекс s опускаем). Имеем

∬ _DN i (x, y) { ^μ _∆ ^b _th n -σ[_∂^∂_x(T b ^∂ _∂ ^h _x ⁿ ) +_∂^∂_y(T b ^∂ _∂ ^h _y ⁿ ) -k b ^hn _m ^- _b ^Hn ] -Q b } dσ+

+∫_sN i (x, y) [σ(T b ^∂hn +β b h n ) -A b ] ds =0,

∬DNi(x, y) {μ∆упtрHn-σ[∂∂x(T∂∂Hxn) +∂∂y(T∂∂Hyn) +kbhnm-bHn-mkii(Hn-Z)] +Q}dσ+ +∫sNi(x,y)[σ(T∂Hn+βHn)-A]ds=0, i= 1,2,…,n, или, после упрощения

∬ _D{N i (x, y) [ ^μ _∆ ^b _th n +σ(k b ^hn _m ^- _b ^Hn ) -Q b ] +σT b ( ^∂ _∂ ^N _x ^i∂ _∂ ^h _x ⁿ + ^∂ _∂ ^N _y ^i∂ _∂ ^h _y ⁿ )} dσ+

∫_S N i (x, y)(σβ b h n - A b )ds = 0,

∬ _D{N i (x,y) [ ^μ _∆ ^уп _t ^р H n -σ(k b ^hn _m ^- _b ^Hn -_m ^kb _b(H n -Z) +Q)] +σT( ^∂ _∂ ^N _x ^i∂ _∂ ^H _x ⁿ + ^∂ _∂ ^N _y ^i∂ _∂ ^H _y ⁿ )} dσ+

∫ Ni(x, y)(σβHn - A)ds = 0, i = 1,2,…, n.

Подставляя вместо функций h_n(x, y, t) и H_n(x, y, t) их разложения (11), получаем системы линейных алгебраических уравнений относительно h j и H j :

Бюллетень науки и практики / Bulletin of Science and Practice Т. 8. №6. 2022

∑_jⁿ₌₁ a ij h j = c i , i = 1,2,…,n,                                        (14)

∑ j ⁿ =1 b ij H j =d i , i= 1,2,…,n,                                   (15)

где aij = ∬D [(μb + σkb) Ni(x, y)Nj(x, y) + σTbq(NiNj)] dσ + σ∫SβbNi(x,y)Nj(x,y)ds, ci = σ ∑jn=1 Hj ∬D kb Ni(x, y)Nj(x, y)dσ + ∬D QbNi(x, y) dσ + ∫S AbNi(x, y) ds, bij = ∬D [(μупр + σkb - σki) Ni(x, y)Nj(x, y) + σTq(NiNj)] dσ + σ∫SβNi(x,y)Nj(x,y)ds, di = σ ∑jn=1 hj ∬D kb Ni(x, y)Nj(x, y)dσ + ∬D (σkn z + Q) Ni(x, y) dσ + ∫S ANi(x, y) ds,

q(N,N)= ^∂Ni∂Nj + ^∂Ni∂Nj . ^j      ∂x ∂x    ∂y ∂y

Система дифференциальных уравнений является нелинейной, так как мощность и водопроводимость потока грунтовых вод зависят от функции h(x, y, t), поэтому для решения этой системы применяется простая итерация. В первом приближении в формуле (4) вместо функции h(x, y, ts) подставим начальные значения этой функции h0(x, y) , а в следующих приближениях - ее значения из предыдущей итерации. В соответствии с этим алгебраические системы (14) и (15) на каждом временном слое решаются до выполнения условий max |h(v) - h(v-1)| ≤ ε, max |H(v) - H(v-1)| ≤ ε max h(ν) -h(ν-1) ≤ ε,     max H (ν) - H (ν-1) ≤ ε, jj           j                               j          j             j где v — номер итерации, ε > 0 — заданное малое число.

Решив задачу (1)-(3) в промежутке [0, T0] при существующих значениях инфильтрации fb(x, y, t) , определяем h(x, y, T0) и H(x, y, T0) . Затем, используя «начальные» условия (8), решаем ретроспективную задачу (7)-(9). Эта задача решается на той же сетке, по тому же алгоритму, что и задача (1)-(3), но в обратном направлении переменной t. В результате получаем поле функции ψ(x,y,t) и по формуле (6) находим управление fb(x, y, t) , а соответствующие УГВ определяются из задачи (1-3). Эта процедура составляет первую итерацию задачи оптимального управления. На последующих итерациях описанная процедура повторяется при значениях управления fb(x, y, t) , полученных из предыдущей итерации. Процесс продолжается до выполнения условия max|h(x , y , T ) - φ(x ,y )| ≤ δ, где δ > 0 — заданное число. i

Выводы

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