Тен (1194485), страница 5
Текст из файла (страница 5)
Было проведено исследование математической модели уровня грунтовых вод в изотропной среде с помощью теории вариационных неравенств – доказана эквивалентность гладкой и вариационной задач, доказано существование и единственность решения вариационной и, как следствие, гладкой задачи. Вариационный метод оказался эффективным методом исследования.
Была проведена численная реализация модели с применением метода конечных элементов. Реализация произведена в двух вариантах: первый на языке C++ без использования сторонних библиотек, и второй – на языке Python с использованием платформы FEniCS. Такое решение обусловлено тем, что на практике в целях численного решения подобных задач с помощью метода конечных элементов не целесообразно создавать реализацию данного метода с самого начала, так как это достаточно сложно и требует времени. Целесообразно использовать готовую библиотеку с реализаций нужных методов. Существует немало таких программных реализаций – эффективных и испытанных инструментов.
Значимость исследования сложно переоценить. В современных условиях рыночной экономики и повсеместного проникновения информационных технологий огромную роль играет математическое моделирование и математические методы, поскольку от них напрямую зависит эффективность той или иной технологии.
СПИСОК ИСПОЛЬЗОВАННЫХ ИСТОЧНИКОВ
1 Михайлов Л.Е. Грунтовые воды. – Л.: изд. ЛПИ. 1982. 40 с.
2 Леонов А.В. Основы гидрогеологии и инженерной геологии: учебное пособие. Томск: Издательство Томского политехнического университета, 2011. 147 с.
3 Бочевер Ф.М., Лапшин Н.Н., Орадовская А.Е. Защита подземных вод от загрязнения. М., Недра, 1979, 256 с.
4 Экланд И., Темам Р. Выпуклый анализ и вариационные проблемы. М.: Мир, 1979. 399 с.
5 Прудников В.Я. Избранные вопросы выпуклого анализа. Хабаровск: Изд. ТОГУ, 2015. 105 с.
6 Я. Бэр, Д. Заславски, С. Ирмей Физико-математические основы фильтрации воды. – М.: Мир, 1971. 452 с.
7 Крукиер Л.А, Шевченко И.В. Сравнение моделей гравитационного режима течения грунтовых вод // Матем. моделирование. 2002. том 14. номер 2. С. 51-60.
8 Кудрявцев Л.Д. Курс математического анализа. В 3 т. – М.: Дрофа, 2003. 351 с.
9 Андреев В.Б. Лекции по методу конечных элементов: Учебное пособие, М.: Издательский отдел факультета ВМиК МГУ им. М.В. Ломоносова, МАКС Пресс, 2010. 2-е изд. 264 с.
10 Зенкевич О. Метод конечных элементов в технике. М.: МИР, 1978. 539 с.
11 Д. Норри, Ж. де Фриз, Ведение в метод конечных элементов: Пер. с англ. – М.: Мир, 1981. 304 с.
12 Л. Сегерлинд, Применение метода конечных элементов: Пер. с англ. – М.: Мир, 1979. 388 с.
13 Susanne C. Brenner, L. Ridgway Scott, The Mathematical Theory of Finite Element Methods, Springer, 2008. 404 p.
14 Власова Е.А., Зарубин В. С., Кувыркин Г. Н. Приближенные методы математической физики. Выпуск XIII, М.: Издательство МГТУ им. Н. Э. Баумана, 2001. 689 с.
15 Абакумов М.В, Гулин А.В. Лекции по численным методам математической физики: учебное пособии. – М.: ИНФРА-М, 2016. 158 с.
16 Масленников В.Н. Дифференциальные уравнения в частных производных: Учебник. – М.: Изд-во РУДН, 1997. 447 с.
17 John M. Stewart, Python for Scientists, O’Railly Learning Series, 2010. 456 p.
ПРИЛОЖЕНИЕ А программный код численной реализации модели
#include <iostream>
#include <stdio.h>
#include <iomanip>
#include <fstream>
#define eps 0.00001
using namespace std;
double getK11(double, double);
double getK12(double, double);
double getK13(double, double);
double getK14(double, double);
double** genKmatrix(double l, double h);
double** sumMatrix(double**a, double **b, int m, int n);
double** T(double**, int n, int m);
int** T(int** mtrx, int m, int n);
double** prod(double**, double**, int m, int n, int q);
double** prod(int**, double**, int m, int n, int q);
double** prod(double**, int**, int m, int n, int q);
double** genMatrix(int, int);
double** genMatrixD(int, int);
double* genVector(int);
double** genMatrixZ(int m, int n);
int** genOmega(int, int, int, int, int);
void show_mtrx(double **, int n);
void show_mtrx(double **, int n, int m);
void show_mtrx(int **, int n, int m);
void show_vec(double*, int);
void show_vec(int*, int);
double** genRK(double, double, double, double*);
double f(double, double); // right f
void payne_irons(double **, double *, int, double, double, double);
double* seidel(int size, double **a, double* b);
bool converge(double *, double *, int);
//
struct Point {
double x;
double y;
};
struct Pp {
int count;
double x;
double y;
};
struct Node {
int num;
int side;
int v1, v2, v3, v4;
Point gv1, gv2, gv3, gv4;
};
double h_ = 0;
int main()
{ //ИСХОДНЫЕ ДАННЫЕ
//////////////////////////////////////////////
int N;
const double A = 1;
double H, h, l;
cout << "enter N: " << endl;
cin >> N;
if (N % 2 != 0) N++;
cout << "N = " << N << endl;
H = A / N;
h = H;
l = H;
h_ = h;
int t = (N + 1)*(N + 1); // кол-во узлов
//////////////////////////////////////////////
//построение "карты" КЭ - конечных элементов - проход по кэ
//массив КЭ-в
Node** map = new Node*[N];
for (int i = 0; i < N; i++)
map[i] = new Node[N];
int count = 1;
// side - 1, если КЭ на границе.
int side = 0;
for (int i = 0; i < N; i++)
for (int j = 0; j < N; j++)
{
// является ли КЭ на границе
if (j == 0 || j == (N - 1) || i == 0 || i == (N - 1) ) side = 1;
else side = 0;
// создание ноды с глобальными коорд и др данными
Point g1 = { i*h, j*h };
Point g2 = { i*h+h, j*h };
Point g3 = { i*h, j*h+h };
Point g4 = { i*h +h , j*h +h };
Node node = { count, side, j + 1 + i*(N + 1), j + 2 + i*(N + 1), j + 1 + (i + 1)*(N + 1), j + 2 + (i + 1)*(N + 1), g1, g2, g3, g4 };
map[i][j] = node;
count++;
}
ofstream goo("nodes.txt");
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++)
{
goo << " #:" << map[i][j].num;
goo << " " << map[i][j].side;
}
goo << endl;
}
goo.close();
// узлы - точки
ofstream goo2("nodes___.txt");
int cnt = 0;
for (int i = 0; i < N + 1; i++) {
for (int j = 0; j < N + 1; j++)
{
if (j == 0 || j == N || i == 0 || i == N) side = 1;
else side = 0;
goo2 << "(" << i <<" , " << j << ")[" <<cnt <<"] " << side << " " ;
cnt++;
}
goo2 << endl;
}
goo2.close();
///////////////////////////////////////////////////////////////////////////////
//сборка глобалбной матрицы жесткости и вектора правой части
double** gK = genMatrixZ(t, t);
double** Km;
int** omega;
int** omegaT;
int s = 0;
//
double** r = genMatrixZ(t, 1);
double** rr;
int c = 0;
cout << "K matrix process starting... " << endl;
for (int i = 0; i < N; i++)
for (int j = 0; j < N; j++)
{
s = map[i][j].side;
switch (s)
{
case 0:// если простой КЭ - не по границе
{
Km = genKmatrix(l, h);
omega = genOmega(N, map[i][j].v1, map[i][j].v2, map[i][j].v3, map[i][j].v4);
omegaT = T(omega, 4, t);
gK = sumMatrix(gK, prod(prod(omegaT, Km, t, 4, 4), omega, t, 4, t), t, t);
c++;
//
double ppp[4];
double ff = f(0, 0);
rr = genRK(1, ff, H, ppp);
r = sumMatrix(r, prod(omegaT, rr, t, 4, 1), t, 1);
//
}break;
case 1://если КЭ по границе
{
Km = genKmatrix(l, h);
omega = genOmega(N, map[i][j].v1, map[i][j].v2, map[i][j].v3, map[i][j].v4);
omegaT = T(omega, 4, t);
gK = sumMatrix(gK, prod(prod(omegaT, Km, t, 4, 4), omega, t, 4, t), t, t);
c++;
//
double ppp[4];
double ff = f(0, 0);
rr = genRK(1, ff, H, ppp);
r = sumMatrix(r, prod(omegaT, rr, t, 4, 1), t, 1);
//
}break;
default:
{
Km = genKmatrix(l, h);
omega = genOmega(N, map[i][j].v1, map[i][j].v2, map[i][j].v3, map[i][j].v4);
omegaT = T(omega, 4, t);
gK = sumMatrix(gK, prod(prod(omegaT, Km, t, 4, 4), omega, t, 4, t), t, t);
c++;
//
double ppp[4];
double ff = f(0, 0);
rr = genRK(1, ff, H, ppp);
r = sumMatrix(r, prod(omegaT, rr, t, 4, 1), t, 1);
//
}break;
}
}
cout << "K matrix is created" << endl;
//вектор столбец в вектор строку - для Зейделя
double* w = new double[t];
for (int i = 0; i < t; i++)
{
w[i] = r[i][0];
}
//
//show_mtrx(gK, t, t);
//show_vec(w, t);
payne_irons(gK, w, N, 1, 1, 1);
double* sol = seidel(t, gK, w);
//double* sol = gauss(gK, w, t);
//show_vec(w, t);
//show_tem(sol, t, N);
//
//вывод матрицы в файл
double dat;
cout << "Seidel process ended" << endl;
ofstream f("matrix.txt");
for (int i = 0; i < t; i++)
{
for (int j = 0; j < t; j++)
{
dat = gK[i][j];
f << setw(10) << dat << " ";
}
f << endl;
}
f.close();
//
//вывод результатов в файл для gnuplot
double data;
int k = 0;
ofstream ff("C:\\coo\\results.txt");
for (int i = 0; i < N + 1; i++)
{
for (int j = 0; j < N + 1; j++)
{
data = sol[k];
k++;
ff << setw(10) << j << "\t" << N - i << "\t" << data << endl;















