Тен (1194485), страница 4

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 4)

– его подпространства. Очевидно, что и не пусто.

Каждая функция может быть однозначно представлена в виде

где – её значение в узле, а – базисная функция (или функция форма, см. рисунок 2) равная единице в узле и нулю во всех других узлах.

Рисунок 3 – Пример функций форм второго порядка

Пусть .

Определение. Назовем приближенным решением [15] задачи (2.2) функцию

такую что

3.5 Вывод расчетных формул

Пусть – произвольный прямоугольный элемент триангуляции T. Пронумеруем его вершины номерами 1, 2, 3 и 4, будем называть эту нумерацию локальной. Нумерация произведена против часовой стрелки. Будем рассматривать приближенное решение безотносительно порядка интерполяционной функции в целях представления формул в более общем виде.

Пусть задается на элементе точками . Обозначим через – функции формы элемента .

Приближённое решение на элементе , представимо в виде

где

– число узлов в конечном элементе,

значение функции в узле,

– функция форма равная единице в узле и нулю в остальных.

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

Представим левую часть равенства (3.1) в матричном виде как

и обозначим

где

Введем обозначения

Теперь интегралы по области можно представить в виде соответствующих сумм интегралов. Имеем

где в свою очередь

Полагая . Используя известную формулу Получим

где

,

аналогично находим

Далее, полагая при , имеем

3.6 Линейный четырехугольный элемент

Интерполяционный полином для прямоугольного элемента с четырьмя узлами имеет вид

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

Рисунок 4 – Четырёхугольный конечный элемент

Пронумерованные узлы и расположение системы координат показаны на рис. 3. В узлах должны быть выполнены следующие условия:

Подстановка этих выражений в формулу (3.2) приводит к системе четырех уравнений, которые могут быть решены относительно :

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

где

Одно из главных различий между этим элементом и симплекс-элементами состоит в том, что градиенты теперь не постоянны, а меняются линейно вдоль одного из координатных направлений. Например,

Градиент в направлении оси постоянен вдоль оси , но меняется линейно по , и, наоборот, постоянна по , но линейно изменяется вдоль оси .

Рисунок 5 – Естественная система координат четырёхугольного конечного элемента

Полученные результаты для прямоугольного элемента могут быть записаны в безразмерной форме с помощью отношений и . Для :

где

Обозначим

тогда функции формы в соотношении (3.3) могут быть представлены в виде произведений безразмерных координат

где

Схематически этот элемент показан на рисунке 5. Введенная только что система координат называется естественной системой координат, потому что координаты при этом изменяются в диапазоне ,

3.7 Вычисление матрицы жёсткости

Как было показано выше

Где вектор функций форм четырёхугольного элемента

Представим это выражение в развернутом виде

аналогично для

С учётом того, что , локальную матрицу жёсткости для элемента можно записать в виде

Теперь перейдем к естественной системе координат.

Так как , и якобиан преобразования имеет вид

и .

Тогда

И, в итоге

Производные функций форм в естественной системе координат имеют следующий вид

Далее вычислим элементы матрицы жесткости.

В силу того, что область интегрирования – это прямоугольник, то

и вычисляя далее аналогично получим локальную матрицу жёсткости

3.8 Вычисление правых частей

Запишем (3.2) в виде

обозначим

– вектор нагрузки конечного элемента .

Вычислим, принимая

Подставляя функции формы и интегрируя получим, окончательно

3.9 Сборка матрицы жёсткости и вектора нагрузок

Сборка глобальной матрицы жёсткости производится с помощь матрицы связности (матрицы связей) [9].

Матрица связности определяется формулой:

где, N – число узлов в сетке.

(здесь знак «=» означает совпадение)

Например, на рис.6 сетка имеет узлов пронумерованных от 1 до 49.

Элемент номер 2 с локальной нумерацией показан на рисунке красным цветом.

Узлы 1, 2, 3, 4 (в локальной нумерации) совпадают с узлами 2, 3, 9 и 10 (в глобальной нумерации) соответственно и не совпадают с остальными.

Тогда матрица имеет вид:

где размер матрицы , (на месте многоточий стоят все нули)

Рисунок 6 – Глобальная и локальная нумерация узлов

Глобальная матрица жёсткости собирается по формуле [9]:

где N – количество конечных элементов.

Правая часть собирается по формуле

3.10 Учет граничных условий

Пусть СЛАУ записано в виде

Граничные условия Дирихле называются главными. Введение ненулевых условий Дирихле осуществляется с помощью процедуры Пейна-Айронса [11]. Если задается условием Дирихле, то диагональный элемент в p-й строке матрицы умножается на очень большое число, а p-й элемент в матрице заменятся произведением того же самого числа на и диагональный элемент.

3.11 Реализация и результаты работы программы

Исследуемая в данной работе модель на практике позволяет вычислить уровень грунтовых вод на некотором участке земли зная этот уровень на её границе. Узнать этот уровень на границе можно проведя измерения с помощью специального оборудования – с некоторым шагом измеряется уровень на границе и полученные данные заносятся в таблицу. Далее можно аппроксимировать полученную табличную функцию или же использовать сразу эти данные в программе для вычисления уровня вод на всем участве.

В данной работе представлены две программные реализации. Первая написана на языке программирования C++ без использования сторонних библиотек, а вторая – на языке Python с использованием библиотеки FEniCS. Программный код численной реализации на C++ приведен в ПРИЛОЖЕНИИ А.

В качестве функции на границе выбрана следующая функция

функция .

На рис. 6, 7 и 8 приведены графики решений полученных в результате выполнения программы на языке C++ при .

Рисунок 7 – График решения при

Рисунок 8 – График решения при

Рисунок 9 – График решения при

Код программы на языке Python

from dolfin import *

# создание триангуляции и опрделение функций форм

mesh = UnitSquareMesh(32, 32)

V = FunctionSpace(mesh, "Lagrange", 1)

# определение границы

def boundary(x):

return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS or x[1] < DOLFIN_EPS or x[1] > 1.0 - DOLFIN_EPS

# определение граничных условий

u0 = Expression("sin(x[0]) + sin(x[1])", degree=2)

bc = DirichletBC(V, u0, boundary)

# определение вариационной задачи

u = TrialFunction(V)

v = TestFunction(V)

f = Constant(1)

a = inner(grad(u), grad(v))*dx

L = f*v*dx

# вычисление

u = Function(V)

solve(a == L, u, bc)

# сохранение в формате vtk

file = File("poisson.pvd")

file << u

# вывод графика решения

plot(u, interactive=True)

Графики решений полученных в результате выполнения программы на языке Python представлены рис. 9, 10 и 11

Рисунок 10 – Результат работы программы Python при

Рисунок 11 – Результат работы программы Python при

Рисунок 12 – Результат работы программы Python при

Для сравнения на рисунках 13, 14 и 15 представлены графики решения задачи полученные в системе Wolfram Mathematica.

Рисунок 13 – График решения, полученный в системе Wolfram Mathematica при

Рисунок 14 – График решения полученный в системе Wolfram Mathematica при

Рисунок 15 – График решения полученный в системе Wolfram Mathematica при

ЗАКЛЮЧЕНИЕ

В результате выполнения выпускной квалификационной работы цель была достигнута в полном объёме и полностью выполнены все поставленные задачи.

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

Была адаптирована известная математическая модель баланса масс подземных вод, поставлена гладкая задача, поставлена вариационная задача.

Характеристики

Список файлов ВКР