Размещено на http://www.allbest.ru/

УДК 532.546

Применение метода конечных элементов с квадратичными базисами к решению задач напорной фильтрации подземных вод

М.У. Мурзакматов, К.А. Исабеков

Аннотация

Рассматриваются два алгоритма приближенного решения уравнения напорной фильтрации, основанные на квадратичных базисных функциях метода конечных элементов (МКЭ) и результаты решения тестовых примеров сравниваются с аналогичными результатами при линейных базисах.

В МКЭ в качестве базисных, наряду с линейными, могут применяться и нелинейные функции [1,2,3]. Практика показывает, что при удачном подборе базисных функций можно существенно повысить точность приближенного решения [4]. напорный фильтрация линейный квадратичный

В данной работе мы рассмотрим применение МКЭ с квадратичными базисами для приближенного решения задачи напорной фильтрации подземных вод:

(1)

(2)

где H=H(x,y) напорная функция; T=T(x,y) водопроводимость пласта; Q=Q(x,y) функция, учитывающая перетоки из ниже и вышележащих пластов; f=f(x,y) функция источников и стоков; =(x,y), =(x,y) заданные функции; /n производная по нормали к границе области; D область фильтрации в плане, Г=D ее граница.

Для приближенного решения задачи (1), (2) область D разбиваем на треугольные элементы и функцию H(x,y) ищем в виде

(3)

Здесь H_i=H(x_i,y_i), N_i(x,y) базисные функции, n количество узлов области D, где определяются значения искомой функции.

Мы рассмотрим два алгоритма, основанные на квадратичных базисных функциях.

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

(4)

где

так называемые l координаты, представляющие отношения площадей элементарных треугольников;

Подставляя в задаче (1), (2) вместо H(x,y) функцию H_n(x,y) из (3), по обобщенному принципу Галеркина получаем систему уравнений

i=1,2,...,n..

После использования формулы Грина в двойном интеграле и разложения (3) функции H_n приходим к следующей системе линейных алгебраических уравнений (СЛАУ)

(5)

где

(6)

Для нахождения коэффициентов и правых частей системы (7) приходится вычислять определенные интегралы от функций, зависящих от lкоординат. Удобнее в качестве переменных интегрирования использовать именно эти координаты. Поэтому необходимо рассмотреть формулы преобразования координат. Декартовы и l координаты связаны соотношениями [1]

x=l_ix_i+ l_jx_j+ l_kx_k ,

y=l_iy_i+ l_jy_j+ l_ky_k ,

1=l_i+ l_j+ l_k ,

где (x_s , y_s) (s=i,j,k) координаты вершин элементов. Исключая из третьего уравнения l_k , т.е. считая l_i и l_j независимыми координатами, имеем

x=(x_i x_k)l_i+(x_jx_k)l_j+x_k ,

y=(y_i y_k)l_i+(y_jy_k)l_j+y_k .

Якобиан преобразования будет равен

,

тогда

При этом производные базисных функций преобразуются по формулам [2]

(7)

а интегралы из (6) по треугольным элементам по формулам

.(8)

Интеграл в правой части (8) вычисляется с помощью квадратурной формулы [1]

где

W_r весовые коэффициенты,

m число точек интегрирования, которое зависит от порядка интерполирования. Порядок интерполирования и весовые коэффициенты определяются суммой показателей степеней координат в каждом члене.

В формулах (6) встречаются три вида подынтегральных функций, поэтому рассмотрим три случая.

Случай 1⁰.

Используя формулы (7), получаем

Случай 2⁰.

Тогда

В этом случае также следует использовать схему второго порядка точности

(12)

Случай 3⁰.

Поэтому

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

(13)

Контурные интегралы в (6) вычисляются по формуле Симпсона.

Алгоритм 2. Суть этого алгоритма состоит в том, что вместо уравнения (1) рассмотрим уравнение с постоянными коэффициентами

и его частные решения

H₁(x)=x²+p₁x+q₁, H₂(y)=y²+p₂y+q₂

используем в качестве базисных функций для решения задачи (1), (2).

Образуем функции формы по формулам

N_i(x,y)=a_i+b_iH₁(x)+c_iH₂(y),

N_j(x,y)=a_j+b_jH₁(x)+c_jH₂(y),

N_k(x,y)=a_k+b_kH₁(x)+c_kH₂(y)

и искомую функцию H(x,y) в области D представим в виде разложения (3). Непосредственной проверкой можно убедиться в том, что функции N_i , N_j , и N_k играют роль функций формы в классическом МКЭ. Применяя принцип Галеркина к задаче (1), (2), получаем СЛАУ (5) с коэффициентами (6). В данном случае, в отличие линейных функций формы и алгоритма 1, к вычислению интегралов (6) привлечение известных квадратурных формул [1] невозможно.

Работа алгоритмов отлажена на решении ряда тестовых примеров. Приведем здесь три из них.

Пример 1. В круге решается задача (1), (2) с линейными и квадратичными базисами при следующих данных: T(x,y)=1, Q(x,y)=0, f(x,y)=4. Искомой функцией является H(x,y)=x²+y²+1. Круг разделен на 6 равных частей, так что элементами являются равносторонние треугольники со стороной, равной, 0.5 (т.е. шаг сетки h=0.5). В табл.1 приведены точные и приближенные значения искомой функции в узлах одного элемента (рис.1), при этом узел i расположен в центре круга, а узлы j и k на окружности. Относительно невысокая точность результатов объясняется величиной шага сетки и тем, что узел q лежит не на самой границе, а на хорде, соединяющей граничные узлы j и k.

Таблица 1 Результаты применения линейных и квадратичных базисов

Пример 2. Область D представляет собой квадрат {1 x 1,1 y 1}, который путем проведения диагоналей последовательно разбивается на 2, 4, 8, 16, 32, 64 элемента. Следовательно, шаг сетки уменьшается от до 0.5. В каждом из шести случаев задача (1), (2) решается МКЭ с линейными и квадратичными базисами при следующих данных: H(x,y)=2(x²+y²)+1, T(x,y)=1, Q(x,y)=0, f(x,y)=8. Результаты счета приведены в табл.2. Поскольку искомая функция и область D обладают центральной и осевой симметрией, в таблице даны результаты, относящиеся к первой четверти координатной плоскости.

Пример 3. Задача (1), (2) решается с данными примера 2 в круге с единичным радиусом (x²+y²1), который последовательно разбивается на 6, 24 и 96 элементов. В первом случае круг разбивается на 6 равных частей. Элементами являются равносторонние треугольники со сторонами, равными 1 (т.е. шаг сетки h=1). Во втором случае проводится окружность радиуса r=0.5. Внутренний круг состоит из 6 равносторонних треугольников (сторона равна 0.5), а круговой пояс разбивается на 18 элементов. В третьем случае проводятся еще две концентрические окружности с радиусами r=0.25 и r=0.75. Внутренний (малый) круг разбивается на 6 треугольников, а три круговые пояса разбиваются соответственно на 18, 30 и 42 треугольников. В табл. 3 приведены результаты, относящиеся к первому квадранту.

В рассмотренных примерах задача (1), (2) решается при постоянных коэффициентах T, Q, f, , и . Это сделано для выяснения влияния базисных функций на точность приближенного решения, так как в таких случаях погрешность осреднения значений указанных функций по элементу отсутствует. В реальных условиях область фильтрации разбивается на крупные фрагменты, в пределах каждого из которых значения гидрогеологических параметров можно считать постоянными. Именно в таких условиях заметно преимущество алгоритма 2 перед линейными базисами, это особенно заметно в примере 1. При уменьшении размеров элементов, как видно из примера 2, результаты применения линейных и квадратичных базисов нивелируются, а интересы дела, наоборот, требуют применения в гидрогеологических расчетах крупных фрагментов.

Квадратичные базисы в алгоритме 2 дают более точные результаты в примере 3. Это объясняется тем, что в этом примере элементами являются равносторонние треугольники, а в примере 2 прямоугольные.

Алгоритм 1 в случае крупной сетки в смысле точности решения занимает среднее положение между другими алгоритмами, а при измельчении сетки теряет преимущество в точности перед линейными базисами (это, повидимому, связано с увеличением погрешности округления в вычислениях). С точки зрения практического применения этот алгоритм явно проигрывает другим, поскольку он требует вдвое больше информации, а размерность матрицы СЛАУ (5) увеличивается четыре раза. Алгоритм 2 требует для своей работы минимальное количество информации, т.е. столько, сколько для линейных базисов и в сочетании с методом фрагментов является наиболее предпочтительным перед другими алгоритмами.

Литература