unit Unit1;
interface
uses
Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
Dialogs, StdCtrls, ExtCtrls;
type
mat2=array[0..100,0..10000] of Real;
TForm1 = class(TForm)
lbl1: TLabel;
lbledt1: TLabeledEdit;
lbledt2: TLabeledEdit;
lbledt3: TLabeledEdit;
lbledt4: TLabeledEdit;
lbledt5: TLabeledEdit;
lbl2: TLabel;
lbl3: TLabel;
lbledt6: TLabeledEdit;
lbledt7: TLabeledEdit;
lbledt8: TLabeledEdit;
lbledt9: TLabeledEdit;
lbl4: TLabel;
lbledt10: TLabeledEdit;
btn1: TButton;
lbledt11: TLabeledEdit;
tmr1: TTimer;
procedure btn1Click(Sender: TObject);
procedure tmr1Timer(Sender: TObject);
private
public
end;
var
Form1: TForm1;
p,v:mat2;
pg,xg,pa, t, dzi, dp, i, pp, v1, a, b, nu, ro, c, l, d, del, tz, vsr, re : Real;
sec, j, k : Integer;
const g=9.81;
implementation
uses Unit2;
{$R *.dfm}
procedure dzit;
begin
if (1-t/tz)<0.05 then
begin
dzi:=1000000;
end
else
begin
if (1-t/tz)<0.95 then
begin
dzi:=Exp(2.3*(7.661175-72.63827*(1-t/tz)+345.7625*exp(2*ln(1-t/tz))-897.8331*exp(3*ln(1-t/tz))+1275.939*exp(4*ln(1-t/tz))-938.8331*exp(5*ln(1-t/tz))+278.8193*exp(6*ln(1-t/tz))))
end
else
begin
dzi:=0.6-0.6*(1-t/tz)
end
end;
end;
function lam(x,m:Integer):Real;
begin
Re := Abs(v[x,m])*d/nu;
if re = 0 then begin Result:=0; end else
if Re <= 2320 then
begin
result := 64 / Re;
end else
begin
result := 0.11*exp(0.25*ln(del/d+68/Re))
end;
end;
function lam1(v:Real):Real;
begin
Re := Abs(v)*d/nu;
if re = 0 then begin Result:=0; end else
if Re <= 2320 then
begin
result := 64 / Re;
end else
begin
result := 0.11*exp(0.25*ln(del/d+68/Re))
end;
end;
function Imin:Real;
begin
result:=p[j-1,k-1]+ro*c*v[j-1,k-1]-0.01*l*lam(j-1,k-1)*ro*v[j-1,k-1]*abs(v[j-1,k-1])/2/d;
end;
function Iplus:Real;
begin
result:=p[j+1,k-1]-ro*c*v[j+1,k-1]+0.01*l*lam(j+1,k-1)*ro*v[j+1,k-1]*abs(v[j+1,k-1])/2/d;
end;
procedure TForm1.btn1Click(Sender: TObject);
begin
pp :=StrToFloat(lbledt1.text)*1000000;
a :=StrToFloat(lbledt2.text);
b :=StrToFloat(lbledt3.text)*3600*3600;
nu :=StrToFloat(lbledt4.text)/1000000;
ro :=StrToFloat(lbledt5.text);
c :=StrToFloat(lbledt9.text);
l :=StrToFloat(lbledt6.text);
d :=StrToFloat(lbledt7.text)/1000;
del:=StrToFloat(lbledt8.text)/1000;
tz :=StrToFloat(lbledt10.text);
pa :=strtofloat(lbledt11.text)*1000000;
vsr:=50;
v1:=Sqrt( (a+(pp-pa)/ro/g) / (0.0827*pi*pi*lam1(vsr)*l/d/16+b*pi*pi*d*d*d*d/16));
for k:=1 to 100 do
begin
vsr:=v1;
v1:=Sqrt( (a+(pp-pa)/ro/g) / (0.0827*pi*pi*lam1(vsr)*l/d/16+b*pi*pi*d*d*d*d/16));
end;
i:=lam1(v1)*v1*v1/2/g/d;
dp:=i*ro*g*l;
for j:= 0 to 100 do
begin
v[j,0]:=v1;
p[j,0]:=dp+pa-j*i*ro*g*l/100;
end;
for k:= 1 to 10000 do
begin
t:=k*0.01*l/c;
dzit;
j:=0;
if Iplus < (pp+ro*g*a) then
begin
v[j,k]:= Sqrt( (64*c*c)/(g*g*b*b*pi*pi*pi*pi*d*d*d*d*d*d*d*d) - 16*(Iplus-pp-ro*g*a)/(ro*g*b*pi*pi*d*d*d*d)) - (8*c)/(g*b*pi*pi*d*d*d*d);
end else begin v[j,k]:=0; end;
p[j,k]:=Iplus + ro*c*v[j,k];
for j:= 1 to 99 do
begin
p[j,k]:=(Imin+Iplus)/2;
v[j,k]:=(Imin-Iplus)/2/ro/c;
end;
j:=100;
if Imin > pa then
begin
v[j,k]:= Sqrt( c*c/dzi/dzi + 2*(Imin-pa)/(ro*dzi) ) - c/dzi;
end
else
begin
v[j,k]:=0;
end;
p[j,k]:= Imin - ro*c*v[j,k];
end;
sec:=0;
Tmr1.Enabled := true;
form2.visible:= true;
end;
procedure TForm1.tmr1Timer(Sender: TObject);
begin
Form2.lnsrsSeries1.Clear;
Form2.lnsrsSeries2.Clear;
for j:=0 to 100 do begin
form2.lnsrsSeries1.AddXY(j*l*0.01,p[j,sec],'',clred);
Form2.lnsrsSeries2.AddXY(j*l*0.01,v[j,sec]*1000000,'',clGreen);
end;
sec:=sec+1;
end;
end.
10.03.2012, 17:00
Неверный результат
Линейный вид
