[KsenIzNN]

Помогите люди добрые! Не могу найти ошибку в программе. Перепробовал уже все что мог - все равно не могу найти.

Суть такова. Рассматривается участок трубопровода. Слева стоит насос, справа закрывается задвижка. Программа строит графики давления красным и скоростей зеленым, решая методом характеристик систему диф.уров в частных производных. Но все формулы перепроверены по многу многу раз. Возможно я как то неправильно воспользовался функциями и где то тут у нас с делфи возникли непонимания. Собственно, "неправильность" проявляется во впадине на графике давления, когда волна отражается.

Очень надеюсь на помощь! Сроки поджимают, от результата зависит очень многое. Заранее спасибо!

Код:

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
    { Private declarations }
  public
    { Public declarations }
  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
   {pg:=p[j,sec-1]/50000;
   xg:=j*7;
   form2.Canvas.Pen.Color:=clBtnFace;
   Form2.Canvas.LineTo(50+round(xg),300-round(pg));
  end;
  form2.Canvas.MoveTo(50,300);
  for j:=0 to 100 do begin
   pg:=p[j,sec]/50000;
   xg:=j*7;
   form2.Canvas.Pen.Color:=clRed;
   form2.Canvas.LineTo(50+round(xg),300-round(pg));
  end;
  form2.lbl1.caption:=FloatToStr((sec-1)*0.01);}
  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.

[KsenIzNN]

добавил файл программы

[icWasya]

>добавил файл программы
А где Unit1.pas Unit1.dfm Unit2.pas Unit2.dfm ?