🚀 go-pugleaf

Witam

 Napisalam w Turbo Pacalu program, którym obliczam wysokosc kolumny w ktorej
suszone sa pewne czasteczki. Liczy wysokosc po ktorej czasteczka bedzie
miala zadana wilgotnosc.

Wszystko jest ok. tylko że w danych tzn drugiej kolumnie tabelki wartosci
nie powinny przekraczac wartosci 0.4 a u mnie dochodza one nawet do 2,033.

Podam treść tego programu jeśli ktos byłby tak uprzejmy i mógł mi to
sprawdzić i wskazać bład byłabym bardzo wdzeczna.

Z góry dziekuje za pomoc i pozdrawiam

Monika


program DMS_suszarka;

uses crt;

const
  row=996;          {gestosc wody}
  ros=1190;         {gestosc suchego materialu}
  g=9.81;           {przyspieszenie ziemskie}
  nmer=9;           {liczba rownan}
  ugx=0.4;          {skladowa osiowa predkosci powietrza}
  wl=0.027;        {strumien materialu7}
  dH0=2500000;      {cieplo parowania wody w 0 st.C}
  cs=4000;          {cieplo wlasciwe ciala stalego}
  cw=4200;          {cieplo wlasciwe wody}
  cpw=1865;         {cieplo wlasciwe pary wodnej nasyconej}
  xkr=0.43;        {wilgotnosc krytyczna materialu}
  p=101300;         {cisnienie atmosferyczne}
  Prd=0.7299;       {liczba Prandla}
  Lambda=0.03442;   {wspolczynnik przewodzenia ciepla czastki}
  D=1.6e-6;         {wspolczynnik dyfuzji}

type
  tabl=array[1..nmer] of extended;

var
   x,                   {poczatek przedzialu calkowania}
   epsilon,             {dokladnosc obliczen}
   h,                   {poczatkowa dlugosc kroku calkowania}
   hmin,                {minimalna dlugosc kroku calkowania}
   h0        :real;     {dlugosc kroku}
   y,                   {wartosc funkcji}
   f         :tabl;    {pochodna funkcji}
   kH,                  {wspolczynnik wnikania ciepla}
   cd,                  {wspolczynnik oporu aerodynamicznego}
   cg,                  {cieplo wlasciwe suchego powietrza}
   cH,                  {cieplo wilgotne}
   csw,                 {cieplo wlasciwe czastki}
   dp,                  {srednica czastki}
   dp0,                 {srednica poczatkowa czastki}
   ff,                  {wspolczynnik korekcyjny}
   ky,                  {wspolczynnik wnikania masy}
   N,                   {liczba czastek}
   Nu,                  {liczba Nusselta}
   pr,                  {preznosc rownowagowa}
   Re,                  {liczba Reynoldsa}
   Le,                  {liczba Lewis'a}
   lg,                  {lepkosc dynamiczna powietrza}
   rog,                 {gestosc powietrza}
   rosw,                {gestosc czastki}
   rosw0,               {gestosc poczatkowa materialu rozpylanego}
   Sp,                  {powierzchnia czastki}
   up,                  {predkosc czastki}
   wd,                  {szybkosc odparowania}
   yr,               {wilgotnosc rownowagowa}
   Yo        : real;    {wilgotnosc poczatkowa}

  i,k        : integer;
  znak       : char;

{zdefiniowanie ukladu rownan rozniczkowych}

procedure func(var x: real; y: tabl; var f: tabl);
  begin
    {dUpx/dh}
      f[1]:=((1-rog/rosw)*g-0.75*cd*up*rog*(y[1]-ugx)/(rosw*dp))/y[1];
    {dUpr/dh}
      f[2]:=(y[2]*y[2]/y[4]-0.75*cd*up*rog*y[2]/(rosw*dp))/y[1];
    {dUpt/dh}
      f[3]:=(-2*y[2]*y[3]/y[4]-0.75*cd*up*rog*y[3]/(rosw*dp))/y[1];
    {dr/dh}
      f[4]:=y[2]/y[1];
    {dY/dh}
      f[5]:=(f[9]*(Yo-y[5]))/y[9]-(wl*f[6])/y[9];
    {dX/dh}
      f[6]:=-6*wd*(1+y[6])/(rosw*dp*y[1]);
    {dTg/dh}

f[7]:=(-(wl/y[9])*((cs+cw*y[6])*f[8]+cw*y[8]*f[6])-(cpw*y[7]+dH0)*f[5])/ch;
    {dTp/dh}

f[8]:=(6*(kH*(y[7]-y[8])-wd*(dH0+cpw*y[7])/(N*Sp)))/(rosw*csw*dp*y[1]);
    {dwg/dh}
      f[9]:=(2*Pi*y[4]*f[4]*ugx*rog)/(1+y[5]);
end;


procedure merson (var x    :real;
                      z    :real;
                  var y,
                      f    :tabl;
                  var h    :real;
                      hmin,
                      e    :real);
{global nmer,func}
var q,r,h3,hs,xs,e5 :real;
             i,j,sw :integer;
              g,s,t :tabl;
     bc,be,bh,br,bx :Boolean;

label 11,22;

begin
if hmin<0 then hmin:=0.01*abs(h);
bh:=true; br:=true; bx:=true; bc:=(e<1); e5:=5*abs(e);
if (((z>x)and(h<=0))or((z<x)and(h>0))) then h:=-h;
11:if bc then
   begin
    xs:=x; for j:=1 to nmer do g[j]:=y[j]
   end;
22:hs:=h; q:=x+h-z; be:=true;
   if ((h>0)and(q>=0))or((h<0)and(q<=0)) then
   begin
    h:=z-x; br:=false
   end;
   h3:=h/3;
for sw:=1 to 5 do
begin
 func(x,y,f);
  for i:=1 to nmer do
  begin
   q:=h3*f[i];
  case sw of
 1:begin r:=q; t[i]:=q; end;
 2:begin r:=0.5*(q+t[i]); end;
 3:begin r:=3*q; s[i]:=r; r:=0.375*(r+t[i]);  end;
 4:begin r:=t[i]+4*q; t[i]:=r; r:=1.5*(r-s[i]); end;
 5:begin r:=0.5*(q+t[i]); q:=abs(2*r-1.5*(q+s[i])); end;
   end{of case};
  y[i]:=g[i]+r; if sw=5 then
  begin
   if bc then
    begin
     r:=abs(y[i]);
    if r<1E-3 then r:=e5 else r:=e5*r;
    if (q>=r) and bx then
      begin
       br:=true; bh:=false; h:=0.5*h;
       if abs(h)<hmin then
        begin
         if h<0 then h:=-hmin else h:=hmin; bx:=false
        end;
       for j:=1 to nmer do y[j]:=g[j];
       x:=xs; goto 22
      end;
    if q>=0.03125*r then be:=false
    end
  end
 end{for i};
     case 5+sw of
  6 :x:=x+h3;
  8 :x:=x+0.5*h3;
  9 :x:=x+0.5*h;
7,10:end{of case};
end{for sw};
if bc then   begin
  if be and bh and br then
    begin
    h:=2*h; bx:=true
    end;
  bh:=true
  end;
  if br then goto 11;
h:=hs
end{of merson};


BEGIN
  for i:=1 to 9 do f[i]:=0;
  clrscr;
  {zdefiniowanie warunkow na wlocie do aparatu}
    y[1]:=0.001;    {predkosc czastki w kierunku osiowym}
    y[2]:=6;        {predkosc czastki w kierunku promieniowym}
    y[3]:=100;      {predkosc czastki w kierunku stycznym}
    y[4]:=0.175;    {odleglosc czastki od osi aparatu}
    y[5]:=0.008;    {wilgotnosc powietrza}
    y[6]:=2.33;     {wilgotnosc materialu}
    y[7]:=523;      {temperatura powietrza}
    y[8]:=293;      {temperatura czastki}
    y[9]:=0.71;     {natezenie strumienia masy powietrza}
    Yo:=y[5];

  {rownania}
  {gestosc czastki}
    rosw0:=((1+y[6])*row)/(y[6]+(row/ros));
  {srednica czastki}
    dp0:=0.0002;
  {liczba czastek}
    N:=6*100/(3600*rosw0*Pi*sqr(dp0)*dp0);
  {ustalenie dlugosci kroku}
    h0:=0.002;
  {ustalenie poczatku przedzialu calkowania}
    x:=0.001;
  {poczatkowa dlugosc kroku calkowania}
    h:=0.1*h0;
  {minimalna dlugosc kroku calkowania}
    hmin:=0.001*h;
  {dokladnosc obliczen}
    epsilon:=1e-6;

    clrscr;

    textcolor(7);

  writeln;
  writeln;
  writeln;
  writeln;
  writeln('                         Projekt suszarki rozpryskowej   ');
  writeln;
  writeln;
  writeln('                            autorzy:');
  writeln;
  writeln('                            -  Maciej Krakowiak');
  writeln('                            -  Damian Wojty˜');
  writeln('                            -  Beniamin Szuszkiewicz');
  writeln('                            -  Dawid Kurc');
  writeln;
  writeln;
  textcolor(4);
  writeln;
  writeln;
  writeln;
  writeln;
  writeln('  Spacja - przewijanie');
  writeln('  Esc - wyj˜cie z programu');
  { delay(60000);  }
          clrscr;
  textcolor(4);
  writeln;
  writeln('                         Projekt suszarki rozpryskowej ');
  writeln;
  textcolor(7);
  writeln('          - model mieszania faz        : stopniowa dyspersja
strugi ');
  writeln('          - typ atomizera              : rozpylacz dyskowy');
  writeln('          - rozpylanie monodyspersyjne');
  writeln('          - wydajnosc                  : 100 kg/h');
  textcolor(4);
  writeln('                               Dane wyjsciowe   ');
  textcolor(7);

  for k:=1 to 80 do write(chr(196));

  textcolor(4);
  write('  Parametry materialu suszonego na wlocie:   ');
  textcolor(7);
  writeln(chr(179));
  writeln('    Strumien dostarczanego materialu         ',chr(179),wl:10:4,'
kg/s');
  writeln('    Gestosc surowki
',chr(179),rosw0:10:0,' kg/m3');
  writeln('    Poczatkowa wilgotnosc materialu
',chr(179),y[6]:10:3,' kg/kg');
  writeln('    Srednica czastki
',chr(179),dp0*1000:10:3,' mm');
  writeln('    Temperatura czastki na wlocie
',chr(179),y[8]-273:10:0,' ',chr(248),'C');
  writeln('    Wilgotnosc krytyczna materialu
',chr(179),xkr:10:3,' kg/kg');
  writeln('    Skladowa osiowa predkosci czastki
',chr(179),y[1]:10:3,' m/s');
  writeln('    Skladowa promieniowa predkosc czastki
',chr(179),y[2]:10:0,' m/s');
  writeln('    Skladowa styczna predkosci czastki
',chr(179),y[3]:10:0,' m/s');
  textcolor(4);
  write('  Parametry powietrza na wlocie:             ');
  textcolor(7);
  writeln(chr(179));
  writeln('    Skladowa osiowa predkosci powietrza
',chr(179),ugx:10:2,' m/s');
  writeln('    Wilgotnosc powietrza na wlocie
',chr(179),y[5]:10:3,' kg/kg');
  writeln('    Temperatura powietrza na wlocie
',chr(179),y[7]-273:10:0,' ',chr(248),'C');
  delay(9000);
  repeat until keypressed;

   clrscr;
   textcolor(4);
        writeln('  h      Upx     Upr     Upt         wg   Promien   Y
X      Tg    Tp ');
        writeln('  m      m/s     m/s     m/s        kg/s     m     kg/kg
kg/kg   řC    řC ');
        for k:=1 to 80 do  write(chr(196));
   textcolor(7);
  i:=1;
  repeat

    {lepkosc powietrza  lg=f(Tg=y[7])}
      lg:=8e-7+7e-8*y[7]-6e-11*sqr(y[7])+3e-14*sqr(y[7])*y[7];
    {gestosc powietrza  lg=f(y[7])}
      rog:=3.3581-0.0121*y[7]+2e-5*sqr(y[7])-1e-8*sqr(y[7])*y[7];
    {cieplo wlasciwe powietrza cg=f(y[7])}
      cg:=960.68+0.2299*y[7]-0.0002*sqr(y[7])+2e-7*sqr(y[7])*y[7];
    {wzgledna predkosc czastek}
      up:=sqrt(sqr(y[1]-ugx)+sqr(y[2])+sqr(y[3]));
    {gestosc czastki}
      rosw:=((1+y[6])*row)/(y[6]+row/ros);
    {srednica czastki}
       if y[6]>=xkr then dp:=dp0*exp((1/3)*ln((rosw0-row)/(rosw-row)));
    {liczba Reynoldsa}
      Re:=dp*rog*up/lg;
    {wspolczynnik oporu aerodynamicznego}
      cd:=(24/Re)+(4/exp((1/3)*ln(Re)));
    {cieplo wilgotne}
      cH:=cg+cpw*y[5];
    {cieplo wlasciwe czastki}
      csw:=cs+cw*y[6];
    {powierzchnia czastki}
      Sp:=Pi*sqr(dp);
    {liczba Nusselta}
      Nu:=2+0.6*sqrt(Re)*exp((1/3)*ln(Prd));
    {wspolczynnik wnikania ciepla}
      kH:=Lambda*Nu/dp;
    {liczba Lewis'a}
      Le:=Lambda/(D*cg*rog);
    {wspolczynnik wnikania masy}
      ky:=kH/(cH*Le);
    {preznosc rownowagowa}
      pr:=exp(23.196-(3816.44)/(y[8]-273.15+227.02));
    {wilgotnosc rownowagowa}
      yr:=0.622*pr/(p-pr);
       if y[6]>=xkr then ff:=1 else ff:=y[6]/xkr;
    {szybkosc odparowania}
      wd:=ky*ff*N*Sp*(yr-y[5]);

    {wywolanie procedury calkowania numerycznego}
      merson(x,x+h0,y,f,h,hmin,epsilon);


writeln(x:5:3,chr(179),y[1]:7:3,y[2]:6,y[3]:7,y[9]:8:4,y[4]:7:3,y[5]:9:5,y[6
]:7:3,y[7]-273:7:1,y[8]-273:6:1);

     inc(i);
     if i=15 then
      begin

           textcolor(4);
           for k:=1 to 80 do  write(chr(196));
           textcolor(7);
          { DELAY(2000);    }
           znak:=readkey;
           if znak=#27 then halt;

        repeat until keypressed ;

   textcolor(7);
        clrscr;
        textcolor(4);
        writeln('  h      Upx     Upr         Upt     wg   Promien   Y
X      Tg    Tp ');
        writeln('  m      m/s     m/s         m/s    kg/s     m     kg/kg
kg/kg   řC    řC ');
        for k:=1 to 80 do  write(chr(196));
        textcolor(7);
        i:=0;
      end;

  {warunek koncowy dzialania programu}
 until (y[6]<=0.065) or (y[7]<=y[8]);


END.