Численное моделирование процессов тепломассопереноса в испаряющемся слое жидкости

          F[nx-1]:=Tm[nx-1,j-1]-kc*Tm[nx,j];    //это свободный член F  для уравнения  нихней строки матрицы   (правая  область или поверхность)

          vc[nx-1]:=0;

        end;

    end;

 

   for i:=2 to nx-1 do      //прямой ход    преобразований к единичной

    begin

    vb[i]:=vb[i]-ka*vc[i-1];

    F[i]:=F[i]-ka*F[i-1];

    va[i]:=0;               //двумя предыдущими действиями  привели va[i] к нулю

 

    vc[i]:=vc[i]/vb[i];

      if i=nx-1 then vc[i]:=0;

    F[i]:=F[i]/vb[i];

    vb[i]:=1;               //двумя предыдущими действиями  привели vb[i] к 1

    end;

 

    for i:=nx-2 downto 1 do      //обратный ход

    begin

    F[i]:=F[i]-vc[i]*F[i+1];

    vc[i]:=0;               // предыдущим действием привели vc[i] к нулю

    end;

 

     for i:=1 to nx-1 do

          Tm[i,j]:=F[i];       //запись результатов в массив температур  в начальный момент (не нулевой, а первый)

 

    //расчет и запись в  массив толщины слоя

    if Tm[nx-1,j-1]>Tm[nx,j] then  //температура  по слою понижается

        S[j]:=S[j-1]+Q*dt-R*dt*(Tm[nx-1,j-1]-Tm[nx,j])/dx

    else S[j]:=S[j-1]+Q*dt-R*dt*(Tm[nx,j]-Tm[nx-1,j-1])/dx;

    if S[j]<0 then  //произошло полное испарение - останов

       begin

          S[j]:=0;

          j:=nt;

          boolAll:=true;

       end;

    progonK.Progon1step(j);

end;

 

procedure Tprogon.Execute;   //температура

var i,j,new_nx:integer;

begin

      Progon1step(1); //для первого  момента времени  методом прогонки

  for j:=2 to nt do     //j-счетчик  по времени

  begin

       Progon1step(j);

      if nx>1 then

      begin

      new_nx:=ceil(xk*S[j] / dx)+1;   //округление - кол-во узлов уместившихся  в слой

      if new_nx<nx then

         begin

            old_nx:=old_nx+(nx-new_nx);    //кол-во узлов не уместившихся в слой

            nx:=new_nx;

            for i:=j to nt do

              begin

                //Tm[nx-1,i]:= Tmk/Tm0;

                Tm[nx,i]:= Tmk/Tm0;

                Kn[nx,i]:= Kn[nx,j];

              end;

            Progon1step(j);            //перерасчет по новому количеству узлов

         end;

         end;

   end;

end;

///////////////////// реализация    Tprogon  конец

 

///////////////////// реализация    TprogonK

constructor TprogonK.Create; //конструктор  потока

var i,j:integer;

begin

      inherited Create;     //конструктор  родителя

      //многие параметры установлены конструктором объкта progon

    setlength(Kn,nx+1,nt+1);  //установка  размерности массивов

      Knk:=Kn0;

      for j:=0 to nt do     //j-счетчик по времени

        begin

            Kn[0,j]:=1;             //изначально концентрация у дна заполняется 1, это равносильно начальной концентрации

            Kn[nx,j]:=Knk/Kn0;      //граничные  условия

        end;

      for i:=0 to nx do     //i-счетчик по координате

        begin Kn[i,0]:=1;         //изначально вся среда (слой) заполняется 1, это равносильно начальной концентрации

        end;

    setlength(va,nx+1);   setlength(vb,nx+1);     setlength(vc,nx+1);

    setlength(F,nx+1);

      Kn[0,1]:=Kn0/Kn0;         Kn[nx,1]:=Knk/Kn0;     //граничные условия температур

      F[0]:= Kn[0,1];      //F[nx]:= Kn[nx,1];       //свободные члены уравнений матрицы

                                           //начальная толщина слоя

      //расчет коэффициентов   уравнений на основе разностной схемы

      kc:=20*(Q-R/dx)*dx/aD;      //остается умножить на разницу  температур

 end;

procedure TprogonK.Progon1step(j:integer);

var i :integer;

begin

      //для j-го момента  времени

       F[0]:= Kn[0,j];      //F[nx]:= Kn[nx,j];       //свободные  члены уравнений матрицы

 

  for i:=1 to nx do     //i-счетчик  по координате   до последнего  узла

    begin

     va[i]:=-1;  vb[i]:=1-kc*(Tm[i,j]-Tm[i-1,j]);   //vc[i]:=kc;        //заполнение  разряженной матрицы

     F[i]:=kc*(Tm[i,j]-Tm[i-1,j])/Kn0;                                  //правые части уравнения

 

      if i=1 then

        begin

          F[1]:= Kn0;       //это свободный член F  для уравнения  верхней строки матрицы     (левая область или дно)

          va[1]:=0;  vc[1]:=0;

          F[1]:= F[1]/vb[1];

          vb[1]:=1;

        end;

        end;

 

   for i:=2 to nx do      //прямой  ход    преобразований к единичной

    begin

    F[i]:=F[i]-va[i]*F[i-1];

    va[i]:=0;               //предыдущими действиями привели va[i] к нулю

    F[i]:=F[i]/vb[i];

    vb[i]:=1;               //предыдущими действиями привели  vb[i] к 1

    end;

    //обратный ход  не  нужен

     for i:=1 to nx do

          Kn[i,j]:=F[i];       //запись результатов в массив  концентраций  в j-момент (не нулевой,  а первый и далее)

end;

///////////////////// реализация    TprogonK  конец

///////////////////    закрытые Переменные  модуля

var   progon: Tprogon;        //    глобальная переменная

//  реализация главной формы

procedure TFormMain.FormCreate(Sender: TObject);

begin           

  Chart1.BottomAxis.Automatic:= false;  Chart1.LeftAxis.Automatic := false;

  Button7Click(Sender);   //очистка

end;

procedure TFormMain.Button7Click(Sender: TObject);   //очистка

begin

     Series1.Clear();

     Series2.Clear();

     Series3.Clear();

     Ypmax:=-High(integer);  //изначально очень маленькие значения

     Ypmin:=High(integer);   //изначально  очень большие значения

     edit2.text:='0';     //время в ноль

      CheckBox1.Checked:=false;

      CheckBox1.Visible:=false;

      boolAll:=false;

end;

 

procedure TFormMain.CheckBox1Click(Sender: TObject);

begin

  if CheckBox1.Checked then OutChartG;

end;

 

procedure TFormMain.Button5Click(Sender: TObject);  //масштаб +

var l0,pr:extended;

begin

pr:=(Chart1.LeftAxis.Maximum - Chart1.LeftAxis.Minimum)/2;   //полуразмах (полуамплитуда)

l0:= Chart1.LeftAxis.Minimum + pr;

Chart1.LeftAxis.Maximum:=l0+pr/2;

Chart1.LeftAxis.Minimum:=l0-pr/2;

end;

procedure TFormMain.Button6Click(Sender: TObject);  //масштаб -

var l0,pr:extended;

begin

pr:=(Chart1.LeftAxis.Maximum - Chart1.LeftAxis.Minimum)/2;   //полуразмах (полуамплитуда)

l0:= Chart1.LeftAxis.Minimum + pr;

Chart1.LeftAxis.Maximum:=l0+pr*2;

Chart1.LeftAxis.Minimum:=l0-pr*2;

end;

 

procedure TFormMain.Button1Click(Sender: TObject);   //   один шаг;  время окончания равно tk+dt

begin

      CheckBox1.Checked:=false;  //снятие галочки

      CheckBox1.Visible:=true;//видно  checkBox

   Series1.Clear();

   Series2.Clear();

   Series3.Clear();  //Диаграмма

 

   xk:=strtofloat(PointOrComma(edit1.text));

   if xk<=0 then xk:=xk_min;                       //защита от ошибки

   tk:=strtofloat(PointOrComma(edit2.text));

 

      nx:=strtoint(PointOrComma(edit3.text));   //  сетка

      nt:=strtoint(PointOrComma(edit4.text));   //  сетка  

 

      Tm0:=strtofloat(PointOrComma(edit10.text));   // =1

      Kn0:=strtofloat(PointOrComma(edit11.text));   // =1

      Tmk:=strtofloat(PointOrComma(edit12.text));   //задаваемая разница температур ( <> 0  )

      Knk:=strtofloat(PointOrComma(edit13.text));   //получаемая разница концентраций растворителя ( <> 0  )

      A:=strtofloat(PointOrComma(edit5.text));

      aD:=strtofloat(PointOrComma(edit6.text));

      Q:=strtofloat(PointOrComma(edit7.text));

      R:=strtofloat(PointOrComma(edit8.text));

     progon:=Tprogon.Create;

     if (tk<=0) or (xk<=0) then

        begin

          progon.Progon1step(1);

          OutChart(1);

          OutChartS(1);

        end

     else

        begin

          progon.Execute;

          OutChart(nt);

          OutChartS(nt);

        end;

 

     //изменение времени на один (10*dt) шаг  если слой еще есть

     if not boolAll then edit2.text:=floattostr(nt*dt+10*dt);

end;

 

procedure TFormMain.OutChart(j:integer);   //отображение на графике температуры  в узлах j-времени

var i:integer;

    x:Extended;

begin

 

  for i:=0 to (nx+old_nx) do

   begin

   x:=i/(nx+old_nx);      //отображаемая  относительная координата x*=x/б0

      if (Ypmax<Tm[i,j]) Then Ypmax:=Tm[i,j];

      if (Ypmin>Tm[i,j]) Then Ypmin:=Tm[i,j];

      if (Ypmax<Kn[i,j]) Then Ypmax:=Kn[i,j];

      if (Ypmin>Kn[i,j]) Then Ypmin:=Kn[i,j];

        if  Ypmin=0 then

              Ypmin:=1;

      if x<S[j]+(dx/2) then   //точка в график вставляется,если  коондината меньше конечной толщины  слоя

        begin

        formMain.Chart1.Series[0].AddXY(x,Tm[i,j],'');

        formMain.Chart1.Series[1].AddXY(x,Kn[i,j],'');

        end;

   end;

   SetChartRange(1,0,Ypmax,Ypmin);

end;

 

procedure TFormMain.OutChartG;   //отображение  на графике положения границы  во времени

var i:integer;

    x:Extended;

begin

     Series1.Clear();

     Series2.Clear();

     Series3.Clear();

     Ypmax:=-High(integer);  //изначально  очень маленькие значения

     Ypmin:=High(integer);   //изначально  очень большие значения

 

  for i:=0 to nt do

   begin

   x:=i/nt;      //отображаемая относительная координата x*=x/б0

 

      if (Ypmax<S[i]) Then Ypmax:=S[i];

      if (Ypmin>S[i]) Then Ypmin:=S[i];

        formMain.Chart1.Series[2].AddXY(x,S[i],''); 

   end;

    SetChartRange(1,0,Ypmax,Ypmin);

end;

 

 procedure TFormMain.OutChartS(j:integer);     //отображение на графике границы слоя    для определенного момента времени

begin  

      formMain.Chart1.Series[2].AddXY(S[j],Ypmax,'');

      formMain.Chart1.Series[2].AddXY(S[j],Ypmin,'');

      SetChartRange(formMain.Chart1.BottomAxis.Maximum,0,Ypmax,Ypmin);

end;

 

procedure TFormMain.EdsChange(Sender: TObject);

begin

Button7Click(Sender);   //очистка, при  изменении исходных данных

end;

 

end.