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

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

полученные, например, с помощью изложенной выше

алгебраической модели.

Численные эксперименты показали, что модель

А.Н.Секундова удобна при моделировании и достаточно хорошо

описывает турбулентные движения несжимаемой среды.

Проинтегрировать уравнение (11.14) можно подобно уравнениям

Навье-Стокса (см. разделы 5, 6).

Еще более совершенной моделью для вычисления поля

турбулентной вязкости является двухпараметрическая K-E

модель. Здесь K – кинетическая энергия турбулентных

пульсаций, а E – скорость диссипации кинетической энергии в

этих пульсациях. Модель задается следующими выражениями:

;

turb







(11.16)

;

eff





































(11.17)

;

eff







































(11.18)

;

EDS

turbK





(11.19)

 

;

EcDcS

turbE





(11.20)

;

turbmoleff





(11.21)

=0,09; с

= 1,44; с

=1,92; σ

=1; σ

=1,3. (11.22)

Начальные и граничные условия обычно принимают в

виде

;;

KcEUcK







(11.23)

где c

, c

– некоторые малые положительные постоянные.

Граничные условия на твердых поверхностях имеют вид

.0 ;0 







ГГ

(11.24)

Здесь n – нормаль к неподвижной твердой поверхности Г.

При использовании выражений (11.16) и (11.20) следует обратить

внимание на то, что обращение в процессе вычислений K и E в

ноль ведет к аварийной остановке счета. Следует предусмотреть

в программе соответствующую защиту.

12. Фрагменты программ

12.1. Система P-U, явная схема, программа на Delphi

В данном подразделе с купюрами приведена программа

решения задачи о распространении жидкости в сосуде. В задаче

рассматривается сосуд с двумя приточными (сверху и снизу) и

одним вытяжным (справа) отверстиями (рис. 3). Задача решается

в системе «давление – скорость».

В листинге 1 приведено частичное описание класса

TForm1 и 3 процедуры.

Первая из них, InitData, инициализирует константы

(шаги по расстоянию h и по времени tau, плотность жидкости

ro, кинематическую вязкость nu, константу c для метода слабой

сжимаемости, погрешность при расчете давления epsP), задает

начальные и граничные условия для полей давления P и

скоростей Ux и Uy, а также задает константы, необходимые для

визуального представления результатов на экране (VisHx и

VisHy определяют шаг сетки в пикселях, а PColKoef

используется для расчета насыщенности цвета при выводе поля

давления). Заметим, что скорости в приточных и вытяжном

отверстиях задаются параболическими профилями; обратите

внимание, как задаются границы приточных и вытяжного

отверстий.

Рис. 3

Вторая процедура, UpdateData, вычисляет значения

полей давления и скоростей в новый момент времени. Для

вычисления нового значения поля давления используется метод

слабой сжимаемости. На твердых стенках в качестве граничного

условия для давления его первую производную принимают

равной 0. А на приточных и вытяжном отверстиях вторую

производную давления принимают равной 0. Далее используется

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

Последняя процедура, Redraw, отвечает за вывод

полученных полей на экран в графическом виде. Поле скоростей

отображается красными стрелками, а поле давлений –

насыщенностью зеленого цвета. Так, слабое давление

отображается темно-зеленым, а сильное – светло-зеленым

цветом.

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

Листинг 1

const

Nx=24; Ny=18;

type

TForm1 = class(TForm)

. . .

h, tau, ro, nu, c, epsP, P0, PColKoef: Real;

Ux, Uy, P, divU: array [1..Nx,1..Ny] of Real;

VisHx, VisHy: Integer;

procedure InitData;

procedure UpdateData;

procedure Redraw;

end;

var

Form1: TForm1;

. . .

procedure TForm1.InitData;

var i, j, i1, i2: Integer;

begin

h:=0.1; tau:=0.0001;

ro:=1; nu:=1; c:=10;

epsP:=0.05; P0:=0;

for i:=1 to Nx do

for j:=1 to Ny do begin

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

P[i,j]:=P0;

end;

i1 := Nx div 4; i2 := Nx div 2 + 1;

for i:=i1+1 to i2-1 do

Uy[i,Ny] := -1.5*(i-i1)*(i2-i);

i1 := Nx div 2; i2 := 3*Nx div 4 + 1;

for i:=i1+1 to i2-1 do

Uy[i,1] := 1.5*(i-i1)*(i2-i);

i1 := Ny div 3; i2 := 2*Ny div 3 + 1;

for i:=i1+1 to i2-1 do

Ux[Nx,i] := 3*(i-i1)*(i2-i);

VisHx:=20; VisHy:=20;

PColKoef:=0.4;

end;

procedure TForm1.UpdateData;

var

NewUx, NewUy: array [1..Nx,1..Ny] of Real;

i, j, i1, i2: Integer;

begin

// Вычисляем P (через div U)

for i:=2 to Nx-1 do

for j:=2 to Ny-1 do

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;

for i:=2 to Nx-1 do

for j:=2 to Ny-1 do

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

{давление на границах}

for i:=2 to Nx-1 do begin

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

P[i,Ny]:=P[i,Ny-1];

end;

for j:=2 to Ny-1 do begin

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

P[Nx,j]:=P[Nx-1,j];

end;

i1 := Nx div 4; i2 := Nx div 2 + 1;

for i:=i1 to i2 do

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

i1 := Nx div 2; i2 := 3*Nx div 4 + 1;

for i:=i1 to i2 do

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

i1 := Ny div 3; i2 := 2*Ny div 3 + 1;

for i:=i1 to i2 do

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

// Вычисляем Ux и Uy

for i:=2 to Nx-1 do

for j:=2 to Ny-1 do

NewUx[i,j] := Ux[i,j] + tau * (

-( (Ux[i,j]+Abs(Ux[i,j]))/2 * (Ux[i,j]-Ux[i-1,j])/h

+ (Ux[i,j]-Abs(Ux[i,j]))/2 * (Ux[i+1,j]-Ux[i,j])/h )

-( (Uy[i,j]+Abs(Uy[i,j]))/2 * (Ux[i,j]-Ux[i,j-1])/h

+ (Uy[i,j]-Abs(Uy[i,j]))/2 * (Ux[i,j+1]-Ux[i,j])/h )

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

+ nu * (Ux[i+1,j]+Ux[i-1,j]+Ux[i,j+1]+Ux[i,j-1]-4*Ux[i,j]) /h/

h );

for i:=2 to Nx-1 do

for j:=2 to Ny-1 do

Ux[i,j]:=NewUx[i,j];

for i:=2 to Nx-1 do

for j:=2 to Ny-1 do

NewUy[i,j] := Uy[i,j] + tau * (

-( (Ux[i,j]+Abs(Ux[i,j]))/2 * (Uy[i,j]-Uy[i-1,j])/h

+ (Ux[i,j]-Abs(Ux[i,j]))/2 * (Uy[i+1,j]-Uy[i,j])/h )

-( (Uy[i,j]+Abs(Uy[i,j]))/2 * (Uy[i,j]-Uy[i,j-1])/h

+ (Uy[i,j]-Abs(Uy[i,j]))/2 * (Uy[i,j+1]-Uy[i,j])/h )

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

+ nu * (Uy[i+1,j]+Uy[i-1,j]+Uy[i,j+1]+Uy[i,j-1]-4*Uy[i,j]) /h/

h );

for i:=2 to Nx-1 do

for j:=2 to Ny-1 do

Uy[i,j]:=NewUy[i,j];

end;

procedure TForm1.Redraw;

var

Bitmap: TBitmap;

hmemdc, hdc1: HDC;

col: TColor;

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

l7, v, c, s: Real;

begin

Bitmap := TBitmap.Create;

hdc1 := GetDC (Handle);

hmemdc := CreateCompatibleDC (hdc1);

Bitmap.Handle := CreateCompatibleBitmap (hdc1, VisHx*(Nx+2),

VisHy*(Ny+2));

with Bitmap.Canvas do begin

Rectangle (VisHx, VisHy, VisHx*Nx+1, VisHy*Ny+1);

for i:=1 to Nx do

for j:=1 to Ny do begin

col := RGB (0, Max (40, Min (240, Round

(120+PColKoef*P[i,j]))), 0);

Brush.Color := col;

Pen.Color := col;

Rectangle (VisHx*i - VisHx div 2, VisHy*j - VisHy div 2,

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

end;

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

for i:=1 to Nx do

for j:=1 to Ny do begin

xa:=VisHx*i;

ya:=VisHy*j;

v := sqrt(sqr(ux[i,j]) + sqr(uy[i,j]));

l7:=0.3*v;

xb := round(xa+ux[i,j]);

yb := round(ya+uy[i,j]);

c := ux[i,j]/(v+0.00001);

s := uy[i,j]/(v+0.00001);

xc := xb-round(l7*(0.865*c+0.5*s));

yc := yb-round(l7*(0.865*s-0.5*c));

xd := xb-round(l7*(0.865*c-0.5*s));

yd := yb-round(l7*(0.865*s+0.5*c));

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

LineTo (xc, yc);

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

end;

SelectObject (hmemdc, Bitmap.Handle);

BitBlt (Canvas.Handle, 0, 0, VisHx*(Nx+2), VisHy*(Ny+2), hmemdc, 0,

0, SRCCOPY);

DeleteDC (hmemdc);

ReleaseDC (Handle, hdc1);

Bitmap.Free;

end;

12.2. Системы P-U и ω-ψ, неявная схема (метод переменных

направлений и расщепления), программа на Delphi

В п.12.2–12.4 решается задача о распространении воздуха

вокруг отдельно стоящего здания. Задача решается в плоскости,

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

(листинг 2) используются две системы для расчета полей

скорости: система «давление – скорость» и система «вихрь –

функция тока». Результат работы программы в системе

«давление – скорость» представлен на рис. 4.

Программа написана на языке Delphi и содержит

несколько классов. Класс TUModel отвечает за расчет поля

скоростей, класс TPModel – за расчет поля давления, класс

TOmPsiModel – за расчет полей функции тока и вихря, а класс

TWholeModel агрегирует перечисленные выше классы.

Большинство переменных в этих классах не нуждаются в

пояснении, комментарии описывают их достаточно подробно.

Поясним назначение матрицы Cells. Эта матрица имеет те же

размеры, что и искомые поля, и каждый элемент этой матрицы

указывает, принадлежит ли соответствующий узел зданию

(kBuild) или нет (kAir).

В конструкторах классов инициализируются переменные

классов, устанавливаются размеры массивов. В конструкторе

класса TWholeModel, кроме того, создаются объекты других

классов, а также вызывается метод InitFields, который

инициализирует рассчитываемые поля. В деструкторе класса

TWholeModel уничтожаются объекты других классов.

Рис. 4

Метод UpdateData класса TPModel рассчитывает поле

давления методом слабой сжимаемости. В этой процедуре

реализована возможность задавать различные граничные условия

на каждой из 4 границ.

Метод UpdatePsi класса TOmPsiModel решает

уравнение Пуассона (8.7) методом верхней релаксации и находит

поле функции тока. Обратите внимание, как вычисляется

константа θ, и как задается поле функции тока на границах.

С помощью процедуры UpdateG класса TWholeModel

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

при небольшой модификации – поля других физических величин,

которые не были рассмотрены в данном пособии. Для этого

достаточно передать процедуре следующие параметры: название

поля (kind: TModelKind), само поле (матрица G),

коэффициент в диффузионной части (alfa: Real; для расчета

полей скорости и вихря, это – кинематическая вязкость), иногда

требуется указать дополнительную информацию (res:

Integer: 0 для поля Ux и 1 для поля Uy). В этой процедуре

реализованы оба метода неявной схемы: метод переменных

направлений и метод расщепления. Параметр meth типа

TMethodKind указывает, какой метод нужно использовать. Эту

процедуру нужно вызывать для каждого этапа в методе

переменных направлений или расщепления, указывая номер

этапа в параметре step. Кроме того, в процедуру можно

передавать различные граничные условия. Код этой процедуры

оптимизирован: с помощью дополнительных переменных

уменьшено число операций деления и умножения, а для

вычисления противоточных производных использовались

условные конструкции.

Метод UpdateData класса TWholeModel выполняет

расчет полей в новый момент времени, вызывая для расчета

каждого поля соответствующий метод соответствующего класса.

Метод Redraw класса TWholeModel выводит поля на

форму. В этой процедуре используется ссылка на объект формы

(Owner: TForm; класс формы называется TMainForm). Это

необходимо для того, чтобы получить доступ к канве формы, и

для того, чтобы определить, какие поля следует выводить. Так,

кнопка TMainForm(Owner).btnP отвечает за вывод поля

давления или вихря, TMainForm(Owner).btnPsi – за вывод

линий тока, TMainForm(Owner).btnU – за вывод поля

скорости. Поля давления и вихря отображаются интенсивностью

зеленого цвета, линии тока – кривыми жирными линиями, а поле

скоростей – стрелками.

Листинг 2

unit UModels;

interface

uses Types, Forms;

type

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

TModelKind = (kUModel, kPModel, kOmModel);

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

TCellKind = (kAir, kBuildings);

// Метод переменных направлений или метод расщепления

TMethodKind = (kVarDirectsMeth=2, kSplitMeth=3);

// Система: давление-скорость или омега-пси

TSystemKind = (kPUSystem, kOmPsiSystem);

TWholeModel = class;

// Описание типа для полей

TArrArrReal = array of array of Real;

PArrArrReal = ^TArrArrReal;

TUModel = class (TObject) // Модель скорости

Nx, Ny: Integer;

h, tau: Real;

kNu: Real; // вязкость

Ux0, Uy0: Real; // скорость ветра

Ux, Uy: TArrArrReal; // поля скоростей по x и по y

Owner: TWholeModel;

constructor Create (AOwner: TWholeModel);

end;

TPModel = class (TObject) // Модель давления

Nx, Ny: Integer;

h, tau: Real;

P0, epsP, kC: Real; // нач.значение P; ошибка P; …

P, deltaP, divU: TArrArrReal; // поле давления P; …

Owner: TWholeModel;

constructor Create (AOwner: TWholeModel);

// Сделать шаг по времени для поля P с граничными условиями Gi

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

end;

TOmPsiModel = class (TObject) // Модель омега-пси

Nx, Ny: Integer;

h, tau: Real;

Om, Psi: TArrArrReal; // поля омега и пси

maxPsi, epsPsi: Real; // макс. Psi; ошибка Psi

Owner: TWholeModel;

constructor Create (AOwner: TWholeModel);

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

procedure UpdatePsi;

end;

PWholeModel = ^TWholeModel;

TWholeModel = class (TObject) // Вся модель

Nx, Ny: Integer; // размер сетки

h, tau, Time: Real; // h, tau и все время