Математика / 5. Математическое моделирование

Новак С.И., Челабчи В.В., Челабчи В.Н.

Одесский национальный морской университет, Украина

Сшивка проекционных решений при моделировании фильтрационных течений

Поле пьезометрического давления при фильтрации жидкости (газа) сквозь пористый или гранулированный материал описывается уравнением (1).

. (1)

Скорость фильтрационного течения определяется согласно Дарси (2).

(2)

где      V_L– скорость фильтрации по направлению L, м/с;

k – проницаемость пористой среды, м²;

m – динамическая вязкость фильтрующейся жидкости (газа), Па·с;

Ф –гидростатическое (пьезометрическое) давления, Па.

Рассматривались граничные условия на поверхности области S: Дирихле (3) и Неймана (условие непротекания) (4).

. (3)

(4)

где N – нормаль к поверхности;

cos(a), cos(b), cos(g) – направляющие косинусы точки на поверхности S.

В соответствии с особенностями методики вся область разбивается на конечное количество локальных областей (ЛО), в которых ищутся решения при условии сшивки решений в локальных областях, которые могут перекрывать друг друга (Рис.1).

С каждой из локальных областей (ЛО), связывается локальная система относительных координат с осями x, y, z чаще всего ориентированными параллельно осям соответственно X, Y, Z.

Назначается несколько типовых (базовых) видов локальных областей со своими локальными координатными системами x_nz, y_nz, z_nz. Форма локальной области, расположение начала локальных координат и размещение узловых точек внутри локальной области в достаточной мере произвольны. Но при сшивке решений в ЛО глобальные координаты точек сшивки решений должны совпадать (Рис.1).

Значения глобальных координат X, Y, Z и значения локальных (относительных) координат x, y, z узловых точек отдельной локальной области (ЛО) связаны зависимостями:

(5)

где - масштабные коэффициенты локальной области с индексом nz,

- глобальные координаты локальной области с индексом nz (глобальные координаты начала локальной системы координат ЛО).

В принципе коэффициенты могут принимать различные значения. Но для удобства выкладок и с целью экономии вычислительных ресурсов имеет смысл принимать одинаковые значения масштабных коэффициентов: где - единый масштабный коэффициент ЛО.

На рис.2 показана схема размещение ЛО в системе глобальных координат.

При переходе от глобальной системе координат к системе локальных координат области получим (6), (7), (8).

, (6)

(7)

. (8)

Функция, аппроксимирующая решение в любой локальной области может иметь различный вид. Главное чтобы она и ее производные были непрерывными и достаточно гладкими. В данной разработке рассматривались аппроксимирующие функции (9).

. (9)

Первые и вторые производные аппроксимирующей функции (2.60) имеют вид (10), (11).

(10)

(11)

При подстановке (10) в уравнение математической модели (2) получим (12). Здесь и далее, для простоты изложения, индекс относящий величину к конкретной точке опущен.

(12)

где

Среднее значение квадрата невязки для узловых точек внутри локальной области (13).

(13)

Для удовлетворения условия (8) в каждой точке (x_msn,y_msn,z_msn ) неподвижной границы S где реализуются условия непротекания должно выполняться (14).

(14)

где

Среднее значение квадрата невязки для точек (x_msn,y_msn,z_msn) на границе S согласно (14) имеет вид(15).

. (15)

Для удовлетворения условия (3) в каждой точке (x_mg,y_mg,z_mg) условной границы G и для реализации условий сшивки должно выполняться (16).

, (16)

где ,

F – значение давления в сходственной точке соседней ЛО.

Среднее значение квадрата невязки для точек на границе G согласно (3) .

, (17)

где mgn – количество узловых точек, где реализуется условие (3) и сшивки.

Для совместного удовлетворения решаемому уравнению, граничным условиям и условиям сшивки решений в смежных ЛО суммируются средние значения квадратов невязок (13), (15), (17) с учетом степени влияния отдельных составляющих (18).

, (18)

где kf – весовые коэффициенты влияния.

Таким образом, выражение для функционала имеет вид:

(19)

В дальнейшем для удобства изложения используется обозначение:

(20)

где R=i+(in+1)×j+(in+1) ×(jn+1) ×k , Rm= in+(in+1)×jn+(in+1) ×(jn+1) ×kn-1.

При 0 £ i £ in, 0 £ j £ jn, 0 £ k £ kn, 0 £ R £ Rm,

Следовательно, выражение (19) с учетом (20) приобретет вид (21).

(21)

Условия минимума для функционала δ_å:

, (22)

где , L=i+(in+1)×j+(in+1) ×(jn+1) ×k.

При 0 £ i £ in, 0 £ j £ jn, 0 £ k £ kn 0 £ L £ Rm.

Таким образом, создается система линейных алгебраических уравнений (СЛАУ) порядка Rm каждое из которых имеет вид:

(23)

Используя обозначения R – индекс столбца матрицы и L – индекс строки матрицы получим выражение (24) для вектора правых частей B(L) и (25) для коэффициентов матрицы A(L,R) системы линейных алгебраических уравнений.

. (24)

(25)

Решение полученных систем линейных алгебраических уравнений осуществляется одним из прямых методов.

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