Филатов Е.Ю., Ясинский Ф.Н. Математическое моделирование течений жидкостей и газов

for jj:=0 to (Ny+1) div ViewStep do begin

i:=ViewStep*ii; j:=ViewStep*jj;

xa := VisHx*i;

ya := VisHy*j;

xb := round(xa+scale*VisHx*Ux[i,j]);

yb := round(ya+scale*VisHy*Uy[i,j]);

v := sqrt(sqr(VisHx*Ux[i,j]) + sqr(VisHy*Uy[i,j]));

l7:=0.4*v;

co := VisHx*Ux[i,j]/(v+0.00001);

si := VisHy*Uy[i,j]/(v+0.00001);

xc := xb-round(scale*l7*(0.865*co+0.5*si));

yc := yb-round(scale*l7*(0.865*si-0.5*co));

xd := xb-round(scale*l7*(0.865*co-0.5*si));

yd := yb-round(scale*l7*(0.865*si+0.5*co));

MoveTo (xa, ya); LineTo (xb, yb);

LineTo (xc, yc);

MoveTo (xb, yb); LineTo (xd, yd);

end;

// Выводим Bitmap на форму; деинициализация

with TMainForm (Owner) do begin

SelectObject (hmemdc, Bitmap.Handle);

BitBlt (Canvas.Handle, x0, y0, Width, Height, hmemdc, 0, 0,

SRCCOPY);

DeleteDC (hmemdc);

ReleaseDC (Handle, hdc1);

Bitmap.Free;

end;

end; // TWholeModel.Redraw

end.

12.3. Система ω-ψ, неявная схема (метод расщепления),

программа на C++

Данная программа (листинг 3) – аналог предыдущей

программы, написанный на C++, за исключением того, что здесь

реализованы лишь система «вихрь – функция тока» и только

метод расщепления, а также изменена структура классов. За

дополнительными разъяснениями обращайтесь к п.12.2.

Результат работы программы представлен на рис. 5.

51

Рис. 5

Листинг 3

#include "windows.h"

#include "math.h"

// Модель: скорости, вихря, тока, температуры, вязкости

enum TModelKind { kUModel, kOmModel, kPsiModel };

// Тип ячейки: воздух, здание

enum TCellKind { kAir, kBuild };

#define nFields 6

#define M_PI 3.14159265358979323846

class CModel // Класс "Модель"

{

public:

int Nx, Ny; // Размеры сетки

double h, tau; // Шаги по расстоянию и по времени

double // Начальные и экстр.значения:

kMu, kLamNu, kNu, kR, // молярная масса, ламинарная вязкость и

вязкость воздуха и унив. газовая пост.

kRo, kP, kT, // плотность, давление и температура воздуха

kUx0, kUy0, // начальная скорость ветра

kEpsPsi; // ошибка Psi

union { // Поля: проекции скорости, вихрь, функция

тока, вспомогательное поле

(double **)Fields[nFields];

struct { double **Ux, **Uy, **Om, **Psi, **P, **NewG; };

};

(char **)Cells; // Карта типов ячеек

double *L, *M; // L и M

52

double F; // Нагрузка на здание

RECT Build; // Координаты здания

// Инициализация

CModel ();

virtual ~CModel();

void AllocFields(); // Выделить память для полей

void FreeFields(); // Освободить память для полей

void Init (int aNx, int aNy, // Инициализировать модель

double ah, double atau, RECT aBuild);

// Расчет шага

void UpdateData(); // Сделать шаг по времени

void FindForce(); // Вычислить силу F

void UpdateP(); // Найти давление P

void UpdatePsi(); // Сделать шаг по времени для поля Psi

// Сделать шаг по времени (часть шага номер step) для модели kind

// (поле G, коэффициент alfa)

void UpdateG (int step, TModelKind kind, double **G, double alfa,

int G1, int G2, int G3, int G4);// (0,0) Г2

// Грагичные условия: ****

// 1 - первого рода; Г1 * * Г3

// 2 - первая производная = 0 ****

// 3 - вторая производная = 0 Г4

};

//////////////////////////////////////////////////////////////////////

// Model Class

// Инициализация

//////////////////////////////////////////////////////////////////////

CModel::CModel()

{

for (int k=0; k<nFields; k++)

Fields[k]=NULL;

Cells=NULL;

L=M=NULL;

} // CModel::CModel

CModel::~CModel()

{

FreeFields();

} // CModel::~CModel

// Выделить память для полей

void CModel::AllocFields()

{

int k, i;

for (k=0; k<nFields; k++)

if (!Fields[k]) {

Fields[k] = new double*[Nx+2];

for (i=0; i<=Nx+1; i++)

Fields[k][i] = new double[Ny+2]; }

53

if (!Cells) {

Cells = new char*[Nx+2];

for (i=0; i<=Nx+1; i++)

Cells[i] = new char[Ny+2]; }

int MN = max(Nx,Ny);

L = new double [MN+2];

M = new double [MN+2];

} // CModel::AllocFields

// Освободить память для полей

void CModel::FreeFields()

{

int k, i;

for (k=0; k<nFields; k++)

if (Fields[k]) {

for (i=0; i<=Nx+1; i++)

delete []Fields[k][i];

delete []Fields[k];

Fields[k]=NULL; }

if (Cells) {

for (i=0; i<=Nx+1; i++)

delete []Cells[i];

delete []Cells;

Cells=NULL; }

if (!L) delete []L;

if (!M) delete []M;

} // CModel::FreeFields

// Инициализировать модель

void CModel::Init (int aNx, int aNy, double ah, double atau, RECT

aBuild)

{

int i, j;

// Инициализируем размеры и константы

FreeFields();

Nx=aNx; Ny=aNy;

h=ah; tau=atau;

Build=aBuild;

AllocFields();

kMu=29*0.001/*кг/моль*/;

kR=8.31/*Дж/моль/К*/;

kLamNu=0.000015;

kNu=1000*kLamNu; /*м2/с*/

// турб вязкость в 1000 раз > ламинарной вязкости 0.015

kP=100000/*Па*/;

kT=295/*К*/;

kRo=kMu*kP/kR/kT/*1,2928/*кг/м3*/;

kUx0=20/*м/с*/; kUy0=0/*м/с*/;

kEpsPsi=0.005;

F=0;

// Инициализируем карту Cells и поля

for (i=0; i<=Nx+1; i++)

54

for (j=0; j<=Ny+1; j++) {

Cells[i][j]=kAir;

Ux[i][j]=0; Uy[i][j]=0;

P[i][j]=0;

Om[i][j]=0;

Psi[i][j] = (j-(Ny+1)/2.0)*h*kUx0; }

for (i=Build.left; i<=Build.right; i++)

for (j=Build.top; j<=Build.bottom; j++) {

Cells[i][j]=kBuild;

Psi[i][j]=0; }

// Инициализируем поля на границах

for (j=0; j<=Ny+1; j++) {

Ux[0][j]=Ux[Nx+1][j]=kUx0;

Uy[0][j]=Uy[Nx+1][j]=kUy0; }

for (i=0; i<=Nx+1; i++) {

Ux[i][0]=Ux[i][Ny+1]=kUx0;

Uy[i][0]=Uy[i][Ny+1]=kUy0; }

} // CModel::Init

//////////////////////////////////////////////////////////////////////

// Расчет шага по времени

//////////////////////////////////////////////////////////////////////

// Сделать шаг по времени

void CModel::UpdateData()

{

int i, j, k;

// Расчет полей

for (i=1; i<=Nx; i++) // Расчет Ux, Uy

for (j=1; j<=Ny; j++)

if (Cells[i][j]==kAir) {

Ux[i][j] = (Psi[i+1][j+1]+Psi[i-1][j+1]-Psi[i+1][j-1]-Psi[i-1]

[j-1])/h/4;

Uy[i][j] =-(Psi[i+1][j+1]+Psi[i+1][j-1]-Psi[i-1][j+1]-Psi[i-1]

[j-1])/h/4; }

for (k=0; k<3; k++) // Расчет Omega

UpdateG (k, kOmModel, Om, kNu, 2,2,3,2);

UpdatePsi(); // Расчет Psi

UpdateP(); // Расчет P для нахождения нагрузки F

FindForce(); // Расчет нагрузки F

} // CModel::UpdateData

// Вычислить силу F

void CModel::FindForce()

{

int i, i1, i2, j, j1, j2;

double FP, FFr, FP1, FP2, tk;

// Вычисление нагрузки через давление

FP=FFr=0;

FP1=FP2=0;

i1=Build.left-1; i2=Build.right+1;

55

for (j=Build.top; j<=Build.bottom; j++) {

if (j==Build.top || j==Build.bottom) tk=0.5; else tk=1;

FP1 += tk*P[i1][j]*h;

FP2 += tk*P[i2][j]*h; }

FP=FP1-FP2;

// Сила трения FFr

FFr=0;

j1=Build.top-1; j2=Build.bottom+1;

for (i=Build.left; i<=Build.right; i++) {

if (i==Build.left || i==Build.right) tk=0.5; else tk=1;

FFr += tk*kRo*kNu*(Ux[i][j1-1]-Ux[i][j1]);

FFr += tk*kRo*kNu*(Ux[i][j2+1]-Ux[i][j2]); }

F=FP+FFr;

} // CModel::FindForce

// Найти давление P

void CModel::UpdateP()

{

int i, j;

double Theta=1.85, EpsP=0.01, newP, e;

// Вычисляем delta P (NewG)

for (i=1; i<=Nx; i++)

for (j=1; j<=Ny; j++)

if (Cells[i][j]==kAir)

NewG[i][j] = -kRo *

(pow(Ux[i+1][j]-Ux[i-1][j],2)/4 + (Ux[i][j+1]-Ux[i][j-1])*

(Uy[i+1][j]-Uy[i-1][j])/2 + pow(Uy[i][j+1]-Uy[i][j-1],2)/4);

// Вычисляем P

do {

e=0;

for (i=1; i<=Nx; i++)

for (j=1; j<=Ny; j++)

if (Cells[i][j]==kAir) {

//newP = (P[i+1][j]+P[i-1][j]+P[i][j+1]+P[i][j-1] - NewG[z]

[i][j]) / 4;

newP = (1-Theta)*P[i][j] + Theta/4 * (P[i+1][j]+P[i-1][j]

+P[i][j+1]+P[i][j-1] - NewG[i][j]);

if (e<fabs(P[i][j]-newP)) e=fabs(P[i][j]-newP);

P[i][j]=newP; }

} while (e>EpsP);

// Вычисляем P на границах

for (j=0; j<=Ny+1; j++) {

P[0][j]=P[1][j];

P[Nx+1][j]=P[Nx][j]; }

for (i=0; i<=Nx+1; i++) {

P[i][0]=kP;

P[i][Ny+1]=kP; } }

for (i=2; i<=Nx-1; i++)

for (j=2; j<=Ny-1; j++)

if (Cells[i][j]==kBuild) {

56

if (Cells[i-1][j]==kAir) P[i][j]=P[i-1][j];

if (Cells[i+1][j]==kAir) P[i][j]=P[i+1][j];

if (Cells[i][j-1]==kAir) P[i][j]=P[i][j-1];

if (Cells[i][j+1]==kAir) P[i][j]=P[i][j+1]; }

} // CModel::UpdateP

// Сделать шаг по времени для поля Psi

void CModel::UpdatePsi()

{

int i, j;

double newPsi, e;

double ksi = pow ( (cos(M_PI/(Nx-1)) + cos(M_PI/(Ny-1))) / 2, 2);

double theta = 2*(1-sqrt(1-ksi))/ksi, theta4=theta/4;//1.85;

// Вычисляем Psi

for (i=1; i<=Nx; i++)

for (j=1; j<=Ny; j++)

if (Cells[i][j]==kAir)

NewG[i][j]=h*h*Om[i][j];

do {

e=0;

for (i=1; i<=Nx; i++)

for (j=1; j<=Ny; j++)

if (Cells[i][j]==kAir) {

newPsi = (1-theta)*Psi[i][j] + theta4 * (Psi[i+1][j]

+Psi[i-1][j]+Psi[i][j+1]+Psi[i][j-1] + NewG[i][j]);

if (e < fabs(Psi[i][j]-newPsi)) e = fabs(Psi[i][j]-newPsi);

Psi[i][j]=newPsi; }

} while (e>kEpsPsi);

// Вычисляем Psi на правой границе

for (j=1; j<=Ny; j++)

Psi[Nx+1][j] = Psi[Nx][j];

} // CModel::UpdatePsi

// Сделать шаг по времени (часть шага номер step) для модели kind

// (поле G, коэффициент alfa)

void CModel::UpdateG (int step, TModelKind kind, double **G,

double alfa, int G1, int G2, int G3, int G4)

{

double tau_hh = tau/h/h, tau_h = tau/h, UTau_h, NuTau_hh;

double aij, bij, cij, fij;

int i, j, i1, j1;

switch (step) {

case 0: // Первая часть шага

if (G1==1) i1=1; else i1=0;

for (j=1; j<=Ny; j++) { // По всем строкам

// Вычисляем прогоночные коэффициенты

// в начале

switch (G1) {

case 1: L[1]=0; M[1]=G[0][j]; break;

case 2: L[1]=1; M[1]=0; break; }

57

// везде

for (i=1; i<=Nx; i++) if (Cells[i][j]==kAir) {

// рядом с другой зоной

if (Cells[i-1][j]==kBuild) {

L[i]=0; M[i]=(Psi[i][j]-Psi[i+1][j])/h/h; }

else if (Cells[i-1][j]!=kAir) {

L[i]=0; M[i]=G[i-1][j]; }

// вычисляем aij, bij, cij, fij

UTau_h=tau_h*Ux[i][j]; NuTau_hh=kNu*tau_hh;

if (Ux[i][j]>0) aij = -UTau_h-NuTau_hh; else aij = -NuTau_hh;

bij = 1 + fabs(UTau_h) + 2*NuTau_hh;

if (Ux[i][j]>0) cij = -NuTau_hh; else cij = UTau_h-NuTau_hh;

fij = G[i][j];

// результат

L[i+1] = -cij/(bij+aij*L[i]);

M[i+1] = (fij-aij*M[i])/(bij+aij*L[i]); }

// Вычисляем новое значение поля G

// в конце

switch (G3) {

case 1: NewG[Nx+1][j] = G[Nx+1][j]; break;

case 2: NewG[Nx+1][j] = -M[Nx+1]/(L[Nx+1]-1); break;

case 3: NewG[Nx+1][j] = -((L[Nx]-2)*M[Nx+1]+M[Nx]) / ((L[Nx]-

2)*L[Nx+1]+1); break; }

// везде

for (i=Nx+1; i>=i1; i--) {

if (i!=Nx+1 && Cells[i][j]==kAir)

NewG[i][j] = L[i+1]*NewG[i+1][j]+M[i+1];

// рядом с другой зоной

if ((Cells[i][j]==kBuild) && (Cells[i-1][j]==kAir))

NewG[i][j] = (Psi[i-1][j]-Psi[i-2][j])/h/h;

else if ((Cells[i][j]!=kAir) && (Cells[i-1][j]==kAir))

NewG[i][j] = G[i][j]; } }

// Сохраняем новое значение поля G

for (i=i1; i<=Nx+1; i++)

for (j=1; j<=Ny; j++)

if (Cells[i][j]==kAir) G[i][j]=NewG[i][j];

break;

case 1: // Вторая часть шага

if (G2==1) j1=1; else j1=0;

for (i=1; i<=Nx; i++) { // По всем столбцам

// Вычисляем прогоночные коэффициенты

// в начале

switch (G2) {

case 1: L[1]=0; M[1]=G[i][0]; break;

case 2: L[1]=1; M[1]=0; break; }

// везде

for (j=1; j<=Ny; j++) if (Cells[i][j]==kAir) {

// рядом с другой зоной

if (Cells[i][j-1]==kBuild) {

L[j]=0; M[j]=(Psi[i][j]-Psi[i][j+1])/h/h; }

else if (Cells[i][j-1]!=kAir) {

L[j]=0; M[j]=G[i][j-1]; }

// вычисляем aij, bij, cij, fij

UTau_h=tau_h*Uy[i][j]; NuTau_hh=kNu*tau_hh;

58

if (Uy[i][j]>0) aij = -UTau_h-NuTau_hh; else aij = -NuTau_hh;

bij = 1 + fabs(UTau_h) + 2*NuTau_hh;

if (Uy[i][j]>0) cij = -NuTau_hh; else cij = UTau_h-NuTau_hh;

fij = G[i][j];

// результат

L[j+1] = -cij/(bij+aij*L[j]);

M[j+1] = (fij-aij*M[j])/(bij+aij*L[j]); }

// Вычисляем новое значение поля G

// в конце

switch (G4) {

case 1: NewG[i][Ny+1] = G[i][Ny+1]; break;

case 2: NewG[i][Ny+1] = -M[Ny+1]/(L[Ny+1]-1); break;

case 3: NewG[i][Ny+1] = -((L[Ny]-2)*M[Ny+1]+M[Ny]) / ((L[Ny]-

2)*L[Ny+1]+1); break; }

// везде

for (j=Ny+1; j>=j1; j--) {

if (j!=Ny+1 && Cells[i][j]==kAir)

NewG[i][j] = L[j+1]*NewG[i][j+1]+M[j+1];

// рядом с другой зоной

if ((Cells[i][j]==kBuild) && (Cells[i][j-1]==kAir))

NewG[i][j] = (Psi[i][j-1]-Psi[i][j-2])/h/h;

else if ((Cells[i][j]!=kAir) && (Cells[i][j-1]==kAir))

NewG[i][j] = G[i][j]; } }

// Сохраняем новое значение поля G

for (i=1; i<=Nx; i++)

for (j=j1; j<=Ny+1; j++)

if (Cells[i][j]==kAir) G[i][j]=NewG[i][j];

// Значения поля G в углах

if ((G1!=1) && (G2!=1)) G[0][0]=(G[0][1]+G[1][0])/2;

if ((G1!=1) && (G4!=1)) G[0][Ny+1]=(G[0][Ny]+G[1][Ny+1])/2;

if ((G3!=1) && (G2!=1)) G[Nx+1][0]=(G[Nx+1][1]+G[Nx][0])/2;

if ((G3!=1) && (G4!=1)) G[Nx+1][Ny+1]=(G[Nx+1][Ny]+G[Nx][Ny+1])/2;

break;

case 2: // Третья часть шага

/* for (i=1; i<=Nx; i++)

for (j=1; j<=Ny; j++)

if (Cells[i][j]==kAir) {

fij=0;

G[i][j] += tau*fij; }*/

break;

}

} // CModel::UpdateG

12.4. Система ω-ψ, неявная схема (метод расщепления),

параллельный вариант, программа на C++

В данном подразделе будет представлен программный

комплекс, состоящий из трех программ, написанный для решения

59

задачи п.12.2, на многопроцессорном компьютере с

использованием технологии распараллеливания.

Этот программный комплекс написан на языке Visual C++

7.0 (как и предыдущая программа) и состоит из трех проектов:

AirAttack, Client и Model. В проекте Model реализуется класс

CModel, который рассчитывает поля скорости, линий тока, вихря

и др. Проект AirAttack – это программа, которая запускается на

главном процессоре (0-й процессор, сервер), который отвечает за

распределение работы между другими процессорами

(клиентами), посылает им данные, получает от них результаты и

выводит результаты на экран. Проект Client – это программа всех

остальных процессоров (клиентов), которая организует расчет

физических полей на своем участке с помощью экземпляра

класса CModel. Распараллеливание реализуется с помощью

интерфейса MPI.

Проект Model – это динамически присоединяемая

библиотека (DLL). Во многом она повторяет предыдущую

программу (п.12.3, листинг 3). Проект Client – это Win32-

приложение, которое, впрочем, почти не использует функций

Windows API. Проект AirAttack – это приложение, написанное с

использованием технологии MFC. Поскольку Visual C++

автоматически создает довольно объемный шаблонный код для

MFC-приложения, ниже (как и в предыдущих листингах)

приводятся лишь выдержки из проекта AirAttack, вставленные

нами в шаблонный код.

Каждый из клиентов использует экземпляр класса CModel

для расчета искомых полей на своем участке. Сервер AirAttack

также использует экземпляр класса CModel. Во-первых, в этой

программе реализована возможность работы в однопроцессорном

режиме. В таком случае всю работу выполняет сервер с помощью

экземпляра класса CModel. Во-вторых, AirAttack использует

экземпляр класса CModel, чтобы сформировать искомые поля на

всей рассчитываемой области из полей на участках, получаемых

им со всех клиентов, которые затем и выводит на экран.

Распараллеливание выполняется следующим образом.

Число клиентов фиксировано (4 процессора; см. рис. 6). Здание

60

Филатов Е.Ю., Ясинский Ф.Н. Математическое моделирование течений жидкостей и газов

Подождите немного. Документ загружается.