Доброго времени суток. Но я себя запутала при задании граничных условий(они на фотогр.) и, соответственно, в задании коэффициентов в прогонке, невязки нарастают (прикрепила скрин). Задача: плоская пластина обтекается ламинарным потоком, на первом начальном участке не должно быть пластины и должно быть условие симметрии, потом поток "натыкается" на эту пластину и образуется погранслой.
Здесь w - тип подзадачи (но одинаково все не считает, что логично): 1 - теплопроводность, 2 - теплопроводность с учетом давления, 3 - система уравнений Навье.
Аппроксимация производных (граничные условия второго рода) конечными разностями и сводится в итоге к равенствам вида u(i;j)=u(i-1;j)
Код:
#include <iostream>
#include <vector>
#include <fstream>
#include "math.h"
using namespace std;
vector<double> progonka(vector<double> b, vector<double> c, vector<double> d, vector<double> f){
int n = f.size();
vector<double> x(n), delta(n), l(n);
delta[0] = -d[0] / c[0];
l[0] = f[0] / c[0];
for (int i = 1; i < n; i++){
delta[i] = -d[i] / (c[i] + b[i] * delta[i - 1]);
l[i] = (f[i] - b[i] * l[i - 1]) / (c[i] + b[i] * delta[i - 1]);
}
x[n - 1] = l[n - 1];
for (int i = n - 2; i >= 0; i--)
x[i] = delta[i] * x[i + 1] + l[i];
return x;
}
int main(){
int XM = 401, YM = 201, XYM, max_time = 4490, w = 3;
XM += 1; YM += 1;
double L = 1.2, H = 0.25, u0 = 1.8, Nu = 9e-4, l = 0.2, Cu = 0.5, eps = 1e-1;
double dx, dy, e, e_first, tau, Ro = 1.225, A = 0.25 / (Ro * u0 * u0);
if (XM >= YM)
XYM = XM;
else
XYM = YM;
vector <double> fui(XM - 1), fuj(YM), fvi(XM), fvj(YM - 1),
tempui(XM - 1), tempuj(YM), tempvi(XM), tempvj(YM - 1), b(XYM), c(XYM), d(XYM);
vector<vector<double> > UN(XM - 1, vector<double>(YM)), U_N(XM - 1, vector<double>(YM)), UN1(XM - 1, vector<double>(YM)),
VN(XM, vector<double>(YM - 1)), V_N(XM, vector<double>(YM - 1)), VN1(XM, vector<double>(YM - 1)),
PN(XM, vector<double>(YM)), P_N(XM, vector<double>(YM)), PN1(XM, vector<double>(YM));
dx = L / (XM - 2);
dy = H / (YM - 2);
tau = Cu * dx / u0;
if (w == 2 || w == 3){
for (int j = 0; j < YM; j++)
for (int i = 0; i < XM - 1; i++)
if (i * dx <= l)
UN[i][j] = u0;
}
for (int q = 1; q <= max_time; q++){
//U
//Прогонка вдоль
for (int j = 1; j < YM - 1; j++){
{
b[0] = 0; c[0] = 1; d[0] = 0; fui[0] = u0;
}
b[XM - 2] = -1; c[XM - 2] = 1; d[XM - 2] = 0; fui[XM - 2] = u0;
for (int i = 1; i < XM - 2; i++){
if (w == 1){
b[i] = tau * (-Nu / (dx * dx));
c[i] = 1 + 2 * tau * Nu / (dx * dx);
d[i] = tau * (-Nu / (dx * dx));
fui[i] = UN[i][j];
}
if (w == 2){
b[i] = tau * (-Nu / (dx * dx) - A * tau / (Ro * dx * dx));
c[i] = 1 + 2 * tau * Nu / (dx * dx) + 2 * A * tau * tau / (Ro * dx * dx);
d[i] = tau * (-Nu / (dx * dx) - A * tau / (Ro * dx * dx));
fui[i] = UN[i][j] - tau / (Ro * dx) * (PN[i + 1][j] - PN[i][j]);
}
if (w == 3){
b[i] = tau * (-Nu / (dx * dx) - A * tau / (Ro * dx * dx) - UN[i][j] / dx);
c[i] = 1 + 2 * tau * Nu / (dx * dx) + 2 * A * tau * tau / (Ro * dx * dx);
d[i] = tau * (-Nu / (dx * dx) - A * tau / (Ro * dx * dx) + UN[i][j] / dx);
fui[i] = UN[i][j] - tau / (Ro * dx) * (PN[i + 1][j] - PN[i][j]);
}
}
tempui = progonka(b, c, d, fui);
for (int i = 0; i < XM - 1; i++)
U_N[i][j] = tempui[i];
}
for (int i = 0; i < XM - 1; i++){
U_N[i][YM - 1] = U_N[i][YM - 2];
if (i * dx <= l){
U_N[i][0] = U_N[i][1];
}
else{
U_N[i][0] = 0;
}
}
// Расчет давления
for (int i = 1; i < XM - 1; i++)
for (int j = 1; j < YM - 1; j++)
P_N[i][j] = PN[i][j] - tau * A * (U_N[i][j] - U_N[i - 1][j]) / dx;
UN1 = U_N;
//U
// Прогонка поперек
for (int i = 1; i < XM - 2; i++){
b[0] = 0; c[0] = -1; d[0] = 1; fuj[0] = 0;
b[YM - 1] = -1; c[YM - 1] = 1; d[YM - 1] = 0; fuj[YM - 1] = 0;
for (int j = 1; j < YM - 1; j++){
if (w == 3){
b[j] = tau * (-Nu / (dy * dy) - 1 / (4 * dy) * (VN[i][j - 1] + VN[i + 1][j - 1]));
c[j] = 1 + 2 * tau * Nu / (dy * dy) + tau / (4 * dy) * (VN[i][j] + VN[i + 1][j] - VN[i][j - 1] - VN[i + 1][j - 1]);
d[j] = tau * (-Nu / (dy * dy) + 1 / (4 * dy) * (VN[i][j] + VN[i + 1][j]));
fuj[j] = U_N[i][j];
} else {
b[j] = tau * (-Nu / (dy * dy));
c[j] = 1 + 2 * tau * Nu / (dy * dy);
d[j] = tau * (-Nu / (dy * dy));
fuj[j] = U_N[i][j];
}
}
tempuj = progonka(b, c, d, fuj);
for (int j = 0; j < YM; j++)
UN1[i][j] = tempuj[j];
}
for (int j = 0; j < YM; j++)
UN1[XM - 2][j] = UN1[XM - 3][j];
//V
// Прогонка вдоль
for (int j = 1; j < YM - 2; j++){
b[0] = 0; c[0] = 1; d[0] = 1; fvi[0] = 0;
b[XM - 1] = -1; c[XM - 1] = 1; d[XM - 1] = 0; fvi[XM - 1] = 0;
for (int i = 1; i < XM - 1; i++){
if (w == 3){
b[i] = tau * (-Nu / (dx * dx) - 1 / (4 * dx) * (UN1[i - 1][j + 1] + UN1[i - 1][j]));
c[i] = 1 + 2 * tau * Nu / (dx * dx) + tau / (4 * dx) * (UN1[i][j + 1] + UN1[i][j] - UN1[i - 1][j + 1] - UN1[i - 1][j]);
d[i] = tau * (-Nu / (dx * dx) + 1 / (4 * dx) * (UN1[i][j + 1] + UN1[i][j]));
fvi[i] = VN[i][j];
} else {
b[i] = tau * (-Nu / (dx * dx));
c[i] = 1 + 2 * tau * Nu / (dx * dx);
d[i] = tau * (-Nu / (dx * dx));
fvi[i] = VN[i][j];
}
}
tempvi = progonka(b, c, d, fvi);
for (int i = 0; i < XM; i++)
V_N[i][j] = tempvi[i];
}
for (int i = 0; i < XM; i++)
V_N[i][YM - 1] = 0;
VN1 = V_N;
//V
// Прогонка поперек
for (int i = 1; i < XM - 1; i++){
b[0] = 0; c[0] = 1; d[0] = 0; fvj[0] = 0;
b[YM - 2] = -1; c[YM - 2] = 1; d[YM - 2] = 0; fvj[YM - 2] = 0;
for (int j = 1; j < YM - 2; j++){
!if (w == 1)....!
if (w == 3){
b[j] = tau * (-Nu / (dy * dy) - A * tau / (Ro * dy * dy) - V_N[i][j] / dy);
c[j] = 1 + 2 * tau * Nu / (dy * dy) + 2 * A * tau * tau / (Ro * dy * dy);
d[j] = tau * (-Nu / (dy * dy) - A * tau / (Ro * dy * dy) + V_N[i][j] / dy);
fvj[j] = V_N[i][j] - tau / (Ro * dy) * (P_N[i][j + 1] - P_N[i][j]);
}
}
tempvj = progonka(b, c, d, fvj);
for (int j = 0; j < YM - 1; j++)
VN1[i][j] = tempvj[j];
}
for (int j = 0; j < YM - 1; j++)
VN1[XM - 1][j] = VN1[XM - 2][j];
e = -1;
// Расчет давления
for (int i = 1; i < XM - 1; i++){
for (int j = 1; j < YM - 1; j++){
PN1[i][j] = P_N[i][j] - tau * A * (VN1[i][j] - VN1[i][j - 1]) / dy;
e = max(fabs(UN[i][j] - UN1[i][j]), e);
}
}
if (q == 1)
e_first = e;
e = e / e_first;
UN = UN1;
VN = VN1;
PN = PN1;
cout << "q = " << q << " eps = " << e << endl;
f1out << q << " " << e << endl;
if (e < eps)
break;
}
return 0;
}