Конечно-разностные схемы моделирования распространения волн
Дипломная работа - Компьютеры, программирование
Другие дипломы по предмету Компьютеры, программирование
еобходимое условие устойчивости явного численного решения некоторых дифференциальных уравнений в частных производных. Как следствие, во многих компьютерных симуляциях временной шаг должен быть меньше определённого значения, иначе результаты будут неправильными[1].
Критерий КФЛ применяется к гиперболическим уравнениям. В двумерном случае он будет иметь вид:
(8)
2.6Условие устойчивости
Параболический характер данной схемы так же отражается и на условии устойчивости. В данном случае оно выгляди как
(10)
где - шаг по времени, h - шаг по пространственным осям, g - ускорение силы тяжести, H - глубина.
Исследует нашу схему на устойчивость.
Уравнение 1.
Заменим .
Для этого должно выполняться условие
Выразим отсюда ?t:
Уравнение 2.
Заменим .
Для этого должно выполняться условие
Выразим отсюда ?t:
Уравнение 3.
Заменим .
Для этого должно выполняться условие
Выразим отсюда ?t:
Условие является самым сильным из полученных. Если усилить его до (при h < 1), то получим искомое неравенство.
2.7Алгоритм решения
1)Задать параметры расчётной сетки.
)Задать начальные значения.
)Рассчитать следующий слой.
)Задать граничные условия на новом слою.
)Вернуться к пункту 3, если слой не последний.
)Вывести полученные данные о высотах и характеристиках в файлы.
)Отобразить поведение волн.
2.8Блок-схема
На рисунке 3 изображена блок-схема алгоритма.
Рисунок 3 - блок-схема алгоритма
2.9Реализация
2.9.1Таблица идентификаторов для программы
В таблице 1 представлены идентификаторы переменных, функций и процедур для языков C/C++, используемых в реализации программ.
Таблица 1 - ИдентификаторыДанныеИдентификатор на языке C/C++Переменные и методы Scheme:Массив значений udouble *** uМассив значений vdouble *** vМассив значений высотdouble *** etaКол-во моментов времениt_numКол-во точек на осяхx_numШаг по времени (секунды)dtШаг по оси (метры)dxГлубина бассейнаHВывод массива значенийoutArray(char *, double ***)Очистка массивовclear()Вывод информации о схемеouput()Переменные и методы SchemeN:Фазовая скоростьcЗадать начальные условияsetInitialConditions()Задать граничные условия на n шагуsetBoundaryConditions(int n)Начать расчётыcalculate()Переменные:Счётчикиi, j, kНачальная координата xx0Начальная координата yy0Начальная высотаeta0aaТестовая схемаtest
2.9.2Реализация на языке С/С++ программы
Представлена реализация алгоритма на языке C/C++.
Реализация явной схемы на языке C/C++
#include
#include
#include Scheme{:*** u;// Массив значений u*** v;// Массив значений v*** eta;// Массив значений etaH;// Глубина бассейнаt_num;// Кол-во моментов времениx_num;// Кол-во точек на оси Xdt;// Шаг по времени (секунды)dx;// Шаг по оси (метры)outArray(char * name, double *** out_array){* fout = fopen(name, "w");(int i = 0; i < t_num; i++){(int j = 0; j < x_num; j++){(int k = 0; k < x_num; k++){(fout, "%f\t", out_array[i][j][k]);
}(fout, "\n");
}(fout, "\n\n");
}(fout);
}:void setInitialConditions() = 0;void setBoundaryConditions(int n) = 0;void calculate() = 0;(int N, int M, double Dx, double H_temp){= Dx;= H_temp;= (dx*dx)/(10*H*10);_num = N;_num = M;= new double**[t_num];= new double**[t_num];= new double**[t_num];(int i = 0; i < t_num; i++){[i] = new double*[x_num];[i] = new double*[x_num];[i] = new double*[x_num];(int j = 0; j < x_num; j++){[i][j] = new double[x_num];[i][j] = new double[x_num];[i][j] = new double[x_num];
}
}
}
~Scheme(){(int i = 0; i < t_num; i++){(int j = 0; j < x_num; j++){[] u[i][j];[] v[i][j];[] eta[i][j];
}[] u[i];[] v[i];[] eta[i];
}[] u;[] v;[] eta;
}clear(){(int i = 0; i < t_num; i++){(int j = 0; j < x_num; j++){(int k = 0; k < x_num; k++){[i][j][k] = 0.0;[i][j][k] = 0.0;[i][j][k] = 0.0;
}
}
}
}ouput(){* foutConfig = fopen("config.txt", "w");(foutConfig, "%d\n%lf\n", x_num, dx);(foutConfig, "%d\n%lf\n", t_num, dt);(foutConfig);("eta.txt", eta);
}
};SchemeN: public Scheme {c;:void setInitialConditions(){x0 = double(x_num)*dx/2.0;y0 = double(x_num)*dx/2.0;eta0 = 1.0;a = 0.2;(int i = 0; i < x_num; i++){(int j = 0; j < x_num; j++){[0][i][j] = 0.0;[0][i][j] = 0.0;[0][i][j] = H + eta0*exp( - a*a*((dx*i-x0)*(dx*i-x0) + (dx*j-y0)*(dx*j-y0)));
}
}
}void setBoundaryConditions(int n){(n t_num - 1);(int i = 1; i < x_num-1; i++){[n][0][i] = eta[n-1][1][i];[n][x_num-1][i] = eta[n-1][x_num-2][i];[n][i][0] = eta[n-1][i][1];[n][i][x_num-1] = eta[n-1][i][x_num-2];[n][0][i] = u[n-1][1][i];[n][x_num-1][i] = u[n-1][x_num-2][i];[n][i][0] = u[n-1][i][1];[n][i][x_num-1] = u[n-1][i][x_num-2];[n][0][i] = v[n-1][1][i];[n][x_num-1][i] = v[n-1][x_num-2][i];[n][i][0] = v[n-1][i][1];[n][i][x_num-1] = v[n-1][i][x_num-2];
}[n][0][0] = eta[n-1][1][1];[n][x_num-1][0] = eta[n-1][x_num-2][1];[n][0][x_num-1] = eta[n-1][1][x_num-2];[n][x_num-1][x_num-1] = eta[n-1][x_num-2][x_num-2];[n][0][0] = u[n-1][1][1];[n][x_num-1][0] = u[n-1][x_num-2][1];[n][0][x_num-1] = u[n-1][1][x_num-2];[n][x_num-1][x_num-1] = u[n-1][x_num-2][x_num-2];[n][0][0] = v[n-1][1][1];[n][x_num-1][0] = v[n-1][x_num-2][1];[n][0][x_num-1] = v[n-1][1][x_num-2];[n][x_num-1][x_num-1] = v[n-1][x_num-2][x_num-2];
}void calculate(){(int i = 0; i < t_num-1; i++){(int j = 1; j < x_num-1; j++){(int k = 1; k < x_num-1; k++){[i+1][j][k] = u[i][j][k] -
((c*c*dt)/2.0)*((eta[i][j+1][k+1] - eta[i][j-1][k+1])/(2.0*dx) +
(eta[i][j+1][k-1] - eta[i][j-1][k-1])/(2.0*dx));[i+1][j][k] = v[i][j][k] -
((c*c*dt)/2.0)*((eta[i][j+1][k+1] - eta[i][j+1][k-1])/(2.0*dx) +
(eta[i][j-1][k+1] - eta[i][j-1][k-1])/(2.0*dx));[i+1][j][k] = eta[i][j][k] +*(-0.5*((u[i][j+1][k+1] - u[i][j-1][k+1])/(2.0*dx) + (u[i][j+1][k-1] -[i][j-1][k-1])/(2.0*dx)) - 0.5*((v[i][j+1][k+1] - v[i][j+1][k-1])/(2.0*dx)
+ (v[i][j-1][k+1] - v[i][j-1][k-1])/(2.0*dx)));(dt*(u[i+1][j][k] + v[i+1][j][k])/dx >= c){("Courant-Friedrichs-Lewy condition failed/n");();
}
}
}(i+1);
}(t_num-1);
}(int N, int M, double Dx, double H_temp): Scheme(N, M, Dx, H_temp), c(sqrt(10*H_temp)) {}
};main(){test(600, 100, 1.0, 1.0);.clear();.setInitialConditions();.calculate();.ouput();("pause");1;
}
2.9.3Исходный код программы для вывода в среде MATLAB
Представлен исходный код программы для вывода в среде MATLAB.
Исходный код программы для вывода в среде MATLAB
clear;config.txt;= config(2);= config(4);= 0:dx:(config(1)-1)*dx;= x;= config(3);= moviein(num);= 1;eta.txt;= 0;= dx*(config(1)-1);r = 1:num= eta(n:n+config(1)-1,:);(x,y,z);([0 xmax 0 xmax -0 2]);('x');('y');('z');(:,r) = getframe;= n+config(1);= 0;= 10;(M, repeat, fps);