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

kMu, kR: Real; // молярная масса и R

kP, kT: Real; // давление и темп. воздуха

kRo: Real; // плотность воздуха

Cells: array of array of TCellKind; // массив типов ячеек

UModel: TUModel; // модель скорости

PModel: TPModel; // модель давления

OmPsiModel: TOmPsiModel; // модель омега-пси

VisHx, VisHy: Integer; // масштаб

Area: TPoint; // координаты области в метрах

bx1, by1, bx2, by2: Real;

Buildings: TRect; // координаты здания в метрах

ViewStep: Integer; // шаг вывода поля скоростей

Enabled, Visible: Boolean; // вкл.ли модель, вывод модели

kMethod: TMethodKind; // метод для неявной схемы 2D

kSystem: TSystemKind; // система

Owner: TForm;

constructor Create (aOwner: TForm; aNx, aNy: Integer;

ah, atau: Real; aArea: TPoint; aBuildings:

TRect; aViewStep: Integer; akMethod: TMethodKind;

akSystem: TSystemKind; aEnabled, aVisible: Boolean);

destructor Destroy; override;

procedure InitFields; // Инициализировать поля

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

// kind (поле G, коэффициент alfa) с использованием метода

// meth. res: 0 для Ux, 1 для Uy

procedure UpdateG (meth: TMethodKind; step: Integer;

kind: TModelKind; G: PArrArrReal;

alfa: Real; res, G1, G2, G3, G4: Integer);

// Сделать шаг по времени (0,0) Г2

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

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

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

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

procedure UpdateData;

procedure Redraw (x0, y0: Integer); // Вывести модель на форму

end;

implementation

uses Windows, Graphics, SysUtils, Math, Unit1;

constructor TUModel.Create (AOwner: TWholeModel);

begin

inherited Create;

// Инициализация полей класса

Owner:=AOwner;

Nx:=Owner.Nx; Ny:=Owner.Ny;

h:=Owner.h; tau:=Owner.tau;

kNu:=15;{в 1000 раз > молекулярной вязкости 0.015}

Ux0:=0{м/с}; Uy0:=0{м/с};

SetLength (Ux, Owner.Nx+2, Owner.Ny+2);

SetLength (Uy, Owner.Nx+2, Owner.Ny+2);

end; // TUModel.Create

41

constructor TPModel.Create (AOwner: TWholeModel);

begin

inherited Create;

// Инициализация полей класса

Owner:=AOwner;

Nx:=Owner.Nx; Ny:=Owner.Ny;

h:=Owner.h; tau:=Owner.tau;

P0:=0{Па};

EpsP:=0.001{Па};

kC:=1000;//400..4600

SetLength (P, Owner.Nx+2, Owner.Ny+2);

SetLength (deltaP, Owner.Nx+2, Owner.Ny+2);

SetLength (divU, Owner.Nx+2, Owner.Ny+2);

end; // TPModel.Create

constructor TOmPsiModel.Create (AOwner: TWholeModel);

begin

inherited Create;

// Инициализация полей класса

Owner:=AOwner;

Nx:=Owner.Nx; Ny:=Owner.Ny;

h:=Owner.h; tau:=Owner.tau;

maxPsi:=0;

epsPsi:=0.05;

SetLength (Om, Owner.Nx+2, Owner.Ny+2);

SetLength (Psi, Owner.Nx+2, Owner.Ny+2);

end; // TOmPsiModel.Create

constructor TWholeModel.Create (aOwner: TForm; aNx, aNy: Integer;

ah, atau: Real; aArea: TPoint; aBuildings: TRect;

aViewStep: Integer; akMethod: TMethodKind;

akSystem: TSystemKind; aEnabled, aVisible: Boolean);

var i, j: Integer;

begin

inherited Create;

// Инициализация полей класса

Owner:=aOwner;

Nx:=aNx; Ny:=aNy;

h:=ah; tau:=atau;

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

kR:=8.31{Дж/моль/К};

kP:=100000{Па};

kT:=295{К};

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

VisHx:=8; VisHy:=8;

Area:=aArea;

Buildings:=aBuildings;

bx1:=aBuildings.Left; by1:=aBuildings.Top;

bx2:=aBuildings.Right; by2:=aBuildings.Bottom;

ViewStep:=aViewStep;

kMethod:=akMethod;

kSystem:=akSystem;

Enabled:=aEnabled;

Visible:=aVisible;

42

// Создаем карту Cells и модели U, P, Om-Psi

SetLength (Cells, Nx+2, Ny+2);

UModel := TUModel.Create (Self);

PModel := TPModel.Create (Self);

OmPsiModel := TOmPsiModel.Create (Self);

InitFields;

end; // TWholeModel.Create

procedure TWholeModel.InitFields;

var i, j: Integer;

begin

Time:=0;

// Инициализируем Ux0, Uy0, P0

UModel.Ux0:=20{м/с}; UModel.Uy0:=0{м/с};

PModel.P0:=0{Па};

if kSystem=kOmPsiSystem then

OmPsiModel.maxPsi := Ux0*Ny*h;

// Инициализируем карту Cells и поля Ux, Uy, P, Om, Psi

for i:=0 to Nx+1 do

for j:=0 to Ny+1 do begin

Cells[i,j]:=kAir;

UModel.Ux[i,j]:=0; UModel.Uy[i,j]:=0;

PModel.P[i,j]:=PModel.P0;

OmPsiModel.Psi[i,j]:=0;

OmPsiModel.Om[i,j]:=0;

end;

for i:=Round(bx1/h) to Round(bx2/h) do

for j:=Round(by1/h) to Round(by2/h) do

Cells[i,j]:=kBuildings;

// Инициализируем поля Ux, Uy, P, Om, Psi на границах

for i:=1 to Ny do begin

UModel.Ux[0,i]:=UModel.Ux0;

UModel.Ux[Nx+1,i]:=UModel.Ux0;

UModel.Uy[0,i]:=UModel.Uy0;

UModel.Uy[Nx+1,i]:=UModel.Uy0;

end;

for j:=0 to Nx+1 do begin

UModel.Ux[j,0]:=UModel.Ux0;

UModel.Ux[j,Ny+1]:=UModel.Ux0;

UModel.Uy[j,0]:=UModel.Uy0;

UModel.Uy[j,Ny+1]:=UModel.Uy0;

end;

end; // TWholeModel.InitFields

destructor TWholeModel.Destroy;

begin

UModel.Free;

PModel.Free;

OmPsiModel.Free;

inherited Destroy;

end; // TWholeModel.Destroy;

procedure TPModel.UpdateData (G1, G2, G3, G4: Integer);

43

const Theta=1.85;

var i, j: Integer; newP, e: Real;

begin

// Вычисляем div U

with Owner.UModel do

for i:=1 to Nx do

for j:=1 to Ny do

if Owner.Cells[i,j]=kAir then

divU[i,j]:=(Ux[i+1,j+1]+Ux[i+1,j-1]-Ux[i-1,j-1]-Ux[i-1,j+1]

+Uy[i-1,j+1]+Uy[i+1,j+1]-Uy[i-1,j-1]-Uy[i+1,j-1])) /4/h;

// Вычисляем P (метод слабой сжимаемости)

for i:=1 to Nx do

for j:=1 to Ny do

if Owner.Cells[i,j]=kAir then

P[i,j] := P[i,j] - kC*tau*divU[i,j];

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

for j:=0 to Ny+1 do

case G1 of

1: P[0,j] := P1;

2: P[0,j] := P[1,j];

3: P[0,j] := 2*P[1,j]-P[2,j];

end;

for i:=1 to Nx do

case G2 of

1: P[i,0] := P2;

2: P[i,0] := P[i,1];

3: P[i,0] := 2*P[i,1]-P[i,2];

end;

for j:=0 to Ny+1 do

case G3 of

1: P[Nx+1,j] := P3;

2: P[Nx+1,j] := P[Nx,j];

3: P[Nx+1,j] := 2*P[Nx,j]-P[Nx-1,j];

end;

for i:=1 to Nx do

case G4 of

1: P[i,Ny+1] := P4;

2: P[i,Ny+1] := P[i,Ny];

3: P[i,Ny+1] := 2*P[i,Ny]-P[i,Ny-1];

end;

with Owner do

for i:=2 to Nx-1 do

for j:=2 to Ny-1 do

if Cells[i,j]=kBuildings then begin

if Cells[i-1,j]=kAir then P[i,j]:=P[i-1,j];

if Cells[i+1,j]=kAir then P[i,j]:=P[i+1,j];

if Cells[i,j-1]=kAir then P[i,j]:=P[i,j-1];

if Cells[i,j+1]=kAir then P[i,j]:=P[i,j+1];

end;

end; // TPModel.UpdateData

procedure TOmPsiModel.UpdatePsi (blOnlyAir: Boolean);

var i, j: Integer; newPsi, Theta, Ksi, e: Real;

44

begin

//Theta:=1.85;

Ksi := Sqr ( (Cos(Pi/(Nx-1)) + Cos(Pi/(Ny-1))) / 2);

Theta := 2 * (1-Sqrt(1-Ksi))/Ksi;

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

repeat

e:=0;

for i:=1 to Nx do

for j:=1 to Ny do

if Owner.Cells[i,j]=kAir then begin

// newPsi := (Psi[i+1,j]+Psi[i-1,j]+Psi[i,j+1]+Psi[i,j-1] +

h*h*Om[i,j]) / 4;

newPsi := (1-Theta)*Psi[i,j]+Theta/4*(Psi[i+1,j]+Psi[i-1,j]

+Psi[i,j+1]+Psi[i,j-1] + h*h*Om[i,j]);

if e<Abs(Psi[i,j]-newPsi) then

e:=Abs(Psi[i,j]-newPsi);

Psi[i,j]:=newPsi;

end;

until e<EpsPsi;

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

for i:=0 to Nx+1 do

Psi[i,0] := -maxPsi/2;

for i:=0 to Nx+1 do

Psi[i,Ny+1] := maxPsi/2;

for j:=1 to Ny do

Psi[0,j] := -maxPsi/2+maxPsi*j/(Ny+1);

// for j:=1 to Ny do

// Psi[Nx+1,j] := Psi[Nx,j];

for j:=1 to Ny do

Psi[Nx+1,j] := 2*Psi[Nx,j]-Psi[Nx-1,j];

end; // TOmPsiModel.UpdatePsi

procedure TWholeModel.UpdateG (meth: TMethodKind; step: Integer; kind:

TModelKind;

G: PArrArrReal; alfa: Real; res, G1, G2, G3, G4: Integer);

var

NewG: TArrArrReal;

L, M: array of Real;

aij, bij, cij, fij, rij, taualfamk_2hh, taumk_2h: Real;

i, j, i1, j1, mk: Integer;

begin

// Устанавливаем размер массивов NewG, Lx, Mx, Ly, My

SetLength (NewG, Nx+2, Ny+2);

SetLength (L, max(Nx,Ny)+2);

SetLength (M, max(Nx,Ny)+2);

if meth = kVarDirectsMeth then mk:=1 else mk:=2;

taualfamk_2hh := tau*alfa*mk/2/h/h;

taumk_2h := tau*mk/2/h;

with UModel do

case step of

1: begin // Шаг 1

45

if G1=1 then i1:=1 else i1:=0;

for j:=1 to Ny do begin // По всем строкам

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

{if kind=kOmModel then with OmPsiModel do

G^[0,j] := -(Psi[0,j]-2*Psi[1,j]+Psi[2,j])/h/h;}

case G1 of

1: begin L[1]:=0; M[1]:=G^[0,j]; end;

2: begin L[1]:=1; M[1]:=0; end;

end;

for i:=1 to Nx do if Cells[i,j]=kAir then begin

if Cells[i-1,j]=kBuildings then begin

L[i]:=0;

if kind=kOmModel then

with OmPsiModel do M[i]:=(Psi[i,j]-Psi[i+1,j])/h/h

else M[i]:=0;

end;

if Ux[i,j]>0 then aij := -taumk_2h*Ux[i,j]-taualfamk_2hh

else aij := -taualfamk_2hh;

bij := 1 + taumk_2h*Abs(Ux[i,j]) + 2*taualfamk_2hh;

if Ux[i,j]>0 then cij := -taualfamk_2hh

else cij := taumk_2h*Ux[i,j]-taualfamk_2hh;

fij := G^[i,j];

if meth = kVarDirectsMeth then begin // Метод перем.напр.

if Uy[i,j]>0 then

fij := fij - tau/2 * Uy[i,j] * (G^[i,j]-G^[i,j-1])/h

+ taualfamk_2hh/mk * (G^[i,j+1]+G^[i,j-1]-2*G^[i,j])

else fij := fij - tau/2 * Uy[i,j]*(G^[i,j+1]-G^[i,j])/h

+ taualfamk_2hh/mk * (G^[i,j+1]+G^[i,j-1]-2*G^[i,j]);

if kind = kUModel then with PModel do

if res = 0 then

fij := fij -tau/2 * (P[i+1,j+1]+P[i+1,j-1]

-P[i-1,j+1]-P[i-1,j-1]) /4/h/kRo

else fij := fij -tau/2 * (P[i+1,j+1]+P[i-1,j+1]

-P[i+1,j-1]-P[i-1,j-1]) /4/h/kRo;

end;

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

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

end;

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

{if kind=kOmModel then with OmPsiModel do

G^[Nx+1,j] := G^[Nx,j];//-(Psi[Nx+1,j]-2*Psi[Nx,j]+Psi[Nx-

1,j])/h/h;}

case G3 of

1: NewG[Nx+1,j]:=G^[Nx+1,j];

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

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

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

end;

i:=Nx;

while i>=i1 do begin

if (Cells[i,j]=kBuildings) and (Cells[i-1,j]=kAir) then

if kind<>kOmModel then NewG[i,j]:=0

else with OmPsiModel do

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

if Cells[i,j]=kAir then

46

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

Dec (i);

end;

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

for i:=i1 to Nx+1 do

for j:=1 to Ny do

if Cells[i,j]=kAir then

G^[i,j]:=NewG[i,j];

end;

2: begin // Шаг 2

if G1=1 then j1:=1 else j1:=0;

for i:=1 to Nx do begin // По всем столбцам

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

{if kind=kOmModel then with OmPsiModel do

G^[i,0] := 2*(Psi[i,1]-Psi[i,0]-UModel.Ux0*h)/h/h;}

case G2 of

1: begin L[1]:=0; M[1]:=G^[i,0]; end;

2: begin L[1]:=1; M[1]:=0; end;

end;

for j:=1 to Ny do if Cells[i,j]=kAir then begin

if Cells[i,j-1]=kBuildings then begin

L[j]:=0;

if kind=kOmModel then

with OmPsiModel do M[j]:=(Psi[i,j]-Psi[i,j+1])/h/h

else M[j]:=0;

end;

if Uy[i,j]>0 then aij := -taumk_2h*Uy[i,j] - taualfamk_2hh

else aij := -taualfamk_2hh;

bij := 1 + taumk_2h*Abs(Uy[i,j]) + 2*taualfamk_2hh;

if Uy[i,j]>0 then cij := -taualfamk_2hh

else cij := taumk_2h*Uy[i,j] - taualfamk_2hh;

fij := G^[i,j];

if meth = kVarDirectsMeth then begin // Метод перем.напр.

if Ux[i,j]>0 then

fij := fij - tau/2 * Ux[i,j] * (G^[i,j]-G^[i-1,j])/h

+ taualfamk_2hh/mk * (G^[i+1,j]+G^[i-1,j]-2*G^[i,j])

else fij := fij - tau/2 * Ux[i,j]*(G^[i+1,j]-G^[i,j])/h

+ taualfamk_2hh/mk * (G^[i+1,j]+G^[i-1,j]-2*G^[i,j]);

if kind = kUModel then with PModel do

if res = 0 then

fij := fij -tau/2 * (P[i+1,j+1]+P[i+1,j-1]

-P[i-1,j+1]-P[i-1,j-1]) /4/h/kRo

else fij := fij -tau/2 * (P[i+1,j+1]+P[i-1,j+1]

-P[i+1,j-1]-P[i-1,j-1]) /4/h/kRo;

end;

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

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

end;

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

{if kind=kOmModel then with OmPsiModel do

G^[i,Ny+1] := 2*(Psi[i,Ny]-Psi[i,Ny+1]+UModel.Ux0*h)/h/h;}

case G4 of

1: NewG[i,Ny+1]:=G^[i,Ny+1];

47

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

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

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

end;

j:=Ny;

while j>=j1 do begin

if (Cells[i,j]=kBuildings) and (Cells[i,j-1]=kAir) then

if kind<>kOmModel then NewG[i,j]:=0

else with OmPsiModel do

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

if Cells[i,j]=kAir then

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

Dec (j);

end;

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

for i:=1 to Nx do

for j:=j1 to Ny+1 do

if Cells[i,j]=kAir then

G^[i,j]:=NewG[i,j];

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

if (G1<>1)and(G2<>1) then G^[0,0]:=(G^[0,1]+G^[1,0])/2;

if (G1<>1)and(G4<>1) then G^[0,Ny+1]:=(G^[0,Ny]+G^[1,Ny+1])/2;

if (G3<>1)and(G2<>1) then G^[Nx+1,0]:=(G^[Nx+1,1]+G^[Nx,0])/2;

if (G3<>1)and(G4<>1) then G^[Nx+1,Ny+1]:=(G^[Nx+1,Ny]

+G^[Nx,Ny+1])/2;

end;

3: if meth = kSplitMeth then begin // Шаг 3 (Метод расщепления)

for i:=1 to Nx do

for j:=1 to Ny do if Cells[i,j]=kAir then begin

if kind = kUModel then with PModel do

if res = 0 then

rij := -(P[i+1,j+1]+P[i+1,j-1]-P[i-1,j+1]-P[i-1,j-

1]) /4/h/kRo

else rij := -(P[i+1,j+1]+P[i-1,j+1]-P[i+1,j-1]-P[i-1,j-

1]) /4/h/kRo

else rij:=0;

G^[i,j] := G^[i,j] + tau*rij;

end;

end; // TWholeModel.UpdateG

procedure TWholeModel.UpdateData;

var steps, i, j, k: Integer;

begin

if not Enabled then Exit;

// Расчет P или Psi {2,2,1,2}

if kSystem=kPUSystem then PModel.UpdateData (1,1,1,1)

else OmPsiModel.UpdatePsi (True);

// Расчет Ux, Uy и Omega

if kMethod=kVarDirectsMeth then steps:=2 else steps:=3;

48

for k:=1 to steps do begin

if kSystem=kPUSystem then begin

UpdateG (kMethod, k,kUModel, @UModel.Ux, UModel.kNu, 0,1,1,1,1);

UpdateG (kMethod, k,kUModel, @UModel.Uy, UModel.kNu, 1,1,1,1,1);

end else begin

UpdateG (kMethod, k, kOmModel, @OmPsiModel.Om, UModel.kNu, 0,

2,2,1,2);

with UModel do with OmPsiModel do

for i:=1 to Nx do

for j:=1 to Ny do

if Cells[i,j]=kAir then begin

Ux[i,j] := (Psi[i,j+1]-Psi[i,j-1])/h/2;

Uy[i,j] := (Psi[i-1,j]-Psi[i+1,j])/h/2;

end;

Time:=Time+tau;

end; // TWholeModel.UpdateData

procedure TWholeModel.Redraw (x0, y0: Integer);

var

Bitmap: TBitmap;

hmemdc, hdc1: HDC;

col: TColor;

xa, ya, xb, yb, xc, yc, xd, yd, i, j, ii, jj: Integer;

scale, fpsi, l7, v, co, si, t1: Real;

rect: TRect;

begin

if not Visible then Exit;

// Подготовка Bitmap'а

with TMainForm (Owner) do begin

Bitmap := TBitmap.Create;

hdc1 := GetDC (Handle);

hmemdc := CreateCompatibleDC (hdc1);

Bitmap.Handle := CreateCompatibleBitmap (hdc1, Width, Height);

end;

with Bitmap.Canvas do begin

col := RGB (0, 0, 0);

Brush.Color := col;

Pen.Color := col;

Rectangle (0, 0, Owner.Width-1, Owner.Height-1);

col := RGB (96, 112, 96);

Brush.Color := col;

Pen.Color := col;

Rectangle (VisHx div 2, VisHy div 2,

VisHx*Nx + VisHx div 2, VisHy*Ny + VisHy div 2);

Rectangle (VisHx div 2 + Round(VisHx*bx1/h),

VisHy div 2 + Round(VisHy*by1/h),

VisHx div 2 + Round(VisHx*bx2/h),

VisHy div 2 + Round(VisHy*by2/h));

49

// Выводим давление P или вихрь Om

if TMainForm(Owner).btnP.Down) then

for i:=1 to Nx do

for j:=1 to Ny do begin

if kSystem=kOmPsiSystem then

t1 := Max (50, Min (170, Round

(120+50*OmPsiModel.Om[i,j])))

else t1 := Max (50, Min (170, Round (120+0.5*(PModel.P[i,j]-

PModel.P0))))

col := RGB (0, Round (t1), 0)

Brush.Color := col;

Pen.Color := col;

Rectangle (VisHx*(i-1) + VisHx div 2,

VisHy*(j-1) + VisHy div 2,

VisHx*i + VisHx div 2, VisHy*j + VisHy div 2);

end

else begin

col := RGB (0, 130, 0);

Brush.Color := col;

Pen.Color := col;

Rectangle (VisHx div 2, VisHy div 2,

VisHx*Nx + VisHx div 2, VisHy*Ny + VisHy div 2);

end;

// Рисуем здание

col := RGB(0,0,0);

Brush.Color := col;

rect.Left := Round(VisHx*bx1/h);

rect.Top := Round(VisHy*by1/h);

rect.Right := Round(VisHx*bx2/h);

rect.Bottom := Round(VisHy*by2/h);

FrameRect (rect);

// Выводим линии тока Psi

if TMainForm(Owner).btnPsi.Down and (kSystem=kOmPsiSystem) then

with OmPsiModel do

for i:=0 to Nx do

for j:=0 to Ny do

for ii:=0 to VisHx-1 do

for jj:=0 to VisHy-1 do

if (Cells[i,j]=kAir) or (Cells[i+1,j+1]=kAir) then

begin

fpsi := ((VisHx-ii)*(VisHy-jj)*Psi[i,j] +

(VisHx-ii)*jj*Psi[i,j+1] +

ii*(VisHy-jj)*Psi[i+1,j] +

ii*jj*Psi[i+1,j+1]) /VisHx/VisHy /maxPsi*20;

if Abs(fpsi-Round(fpsi))<0.2/VisHy*h then

Pixels[VisHx*i+ii,VisHy*j+jj] := RGB(127,255,255);

end;

// Выводим скорость U

if TMainForm(Owner).btnU.Down then with UModel do begin

scale:=0.03*ViewStep;

Pen.Color := RGB (255, 127, 127);

for ii:=0 to (Nx+1) div ViewStep do

50

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

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