Форум программистов
 

Восстановите пароль или Зарегистрируйтесь на форуме, о проблемах и с заказом рекламы пишите сюда - alarforum@yandex.ru, проверяйте папку спам!

Вернуться   Форум программистов > IT форум > Помощь студентам
Регистрация

Восстановить пароль
Повторная активизация e-mail

Купить рекламу на форуме - 42 тыс руб за месяц

Ответ
 
Опции темы Поиск в этой теме
Старый 19.11.2017, 00:46   #1
Odyle
Новичок
Джуниор
 
Регистрация: 19.11.2017
Сообщений: 2
По умолчанию Корректность задания граничных условий (C++)

Доброго времени суток. Но я себя запутала при задании граничных условий(они на фотогр.) и, соответственно, в задании коэффициентов в прогонке, невязки нарастают (прикрепила скрин). Задача: плоская пластина обтекается ламинарным потоком, на первом начальном участке не должно быть пластины и должно быть условие симметрии, потом поток "натыкается" на эту пластину и образуется погранслой.

Здесь 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;
}
Изображения
Тип файла: jpg photo_2017-11-19_00-24-36.jpg (51.1 Кб, 123 просмотров)
Тип файла: jpg photo_2017-11-19_00-35-01.jpg (46.6 Кб, 128 просмотров)
Odyle вне форума Ответить с цитированием
Старый 19.11.2017, 13:17   #2
СтудПом
Форумчанин
 
Регистрация: 08.11.2017
Сообщений: 347
По умолчанию Помощь по математике

Это вам на форум помощи по математике. Программирования не много и задача не сложная, но без понимания анализа и вообще вашей задачи, никто не поможет. Разбираться в коде нет смысла.
СтудПом вне форума Ответить с цитированием
Старый 19.11.2017, 14:20   #3
Odyle
Новичок
Джуниор
 
Регистрация: 19.11.2017
Сообщений: 2
По умолчанию

Я уже разобралась, не актуально
Odyle вне форума Ответить с цитированием
Ответ


Купить рекламу на форуме - 42 тыс руб за месяц

Опции темы Поиск в этой теме
Поиск в этой теме:

Расширенный поиск


Похожие темы
Тема Автор Раздел Ответов Последнее сообщение
Проверка на корректность ввода Arcasha Общие вопросы C/C++ 3 17.10.2014 16:00
Корректность xml данных Utkin Общие вопросы по программированию, компьютерный форум 4 06.07.2013 19:49
Корректность ввода на СИ glebast Помощь студентам 4 14.09.2012 22:42
Корректность скобок! Sport Помощь студентам 3 22.03.2012 20:33
Корректность закачки _Den_1984 Работа с сетью в Delphi 0 15.02.2011 13:30