Метод конечных объемов

[Pphantom]

larts в сообщении #861267 писал(а):

Подскажите, пожалуйста, как это можно исправить.

Может, просто код выложите? Так догадаться о том, что происходит, боюсь, нереально...

[larts]

Пока готовил код для посторонних глаз, нашел косяк, исправил, теперь еще другая проблема. Решение расходится при размере ячейки меньше 1. Правда, этого можно избежать, если умножить коэффициенты в матрице на какое-то число, для шага 0,5, например, на 2, для 0,025 - на 600. (Этот момент там, в самом низу)
Вот код, немного большой, но я старался так сгруппировать, чтоб понятно было..

class vol
{
public:
   unsigned char state;   //byte 0 for w; 1 - e; 2 - n; 3 - s; байт, в котором хранится информация о том, какая сторона ячейки - граничная

   coord a1;      //координаты вершин
   coord a2;
   coord a3;
   coord a4;
   coord c;      //координаты центра масс
   float U;
   float V;
   float Vt;
   float Ut;
   float P;
   int i;
   int j;
   vol(float x1, float y1,float x2, float y2,float x3, float y3,float x4, float y4,int i1,int j1)
   {
      a1.x = x1;
      a1.y = y1;
      a2.x = x2;
      a2.y = y2;
      a3.x = x3;
      a3.y = y3;
      a4.x = x4;
      a4.y = y4;
      i = i1;
      j = j1;
      U = 0;                           
      V = 0;
      Ut = 0;
      Vt = 0;
      P = 1;
      state = 0;
   }
   void CompTemp();            //скорость без учета давления
   void AddBToMatr();            //составление вектора правых частей для СЛАУ
   void Correct();   
   void CompC()
   {   
      float s1 = ((a1.x-a2.x)*(a4.y-a2.y)-(a4.x-a2.x)*(a1.y-a2.y))/2;
      float s2 = -0.5*((a1.x-a3.x)*(a4.y-a3.y)-(a4.x-a3.x)*(a1.y-a3.y));
      float v = 0.5*((a4.x-a1.x)*(a2.y-a3.y)-(a2.x-a3.x)*(a4.y-a1.y));
      c.x = ((a1.x+a2.x+a4.x)/3*s1+(a1.x+a3.x+a4.x)/3*s2)/v;      //сначала ищутся центры масс треугольников, потом 
      c.y = ((a1.y+a2.y+a4.y)/3*s1+(a1.y+a3.y+a4.y)/3*s2)/v;      //по формуле - ц.м. четырехугольника
   }

};
vector<vector<vol>> vec;
void vol::Correct()
{
   float Pw = 0;
   float Pe = 0;
   float Pn = 0;
   float Ps = 0;
   P = x.items[(i)*(ny)+(j)];      //х - это вектор с давлением, который получаем в результате решения СЛАУ
   if((bool)(state&1))
      Pw = P;
   else
      Pw = x.items[(i-1)*(ny)+(j)];
   
   if((bool)(state&2))
      Pe = P; 
   else
      Pe = x.items[(i+1)*(ny)+(j)];         
   if((bool)(state&4))
      Pn = P;
   else
      Pn = x.items[(i)*(ny)+(1+j)];
   if((bool)(state&8))
      Ps = P;  
   else
      Ps = x.items[(i)*(ny)+(j-1)];
   float swy = (a2.x-a1.x);         //компоненты площади
   float ssy = -(a3.x-a1.x);
   float sny = (a4.x-a2.x);
   float sey = -(a4.x-a3.x);

   float swx = -(a2.y-a1.y);
   float ssx = (a3.y-a1.y);
   float snx = -(a4.y-a2.y);
   float sex = (a4.y-a3.y);
   
   float v = 0.5*((a4.x-a1.x)*(a2.y-a3.y)-(a2.x-a3.x)*(a4.y-a1.y));

   if(i == 1)         //на входе  и выходе градиент давления сохраняется
      Pw = 2*P-Pe;
   if(i == nx-2)
      Pe = 2*P-Pw;

   float pe = (Pe+P)/2;
   float pw = (Pw+P)/2;
   float ps = (Ps+P)/2;
   float pn = (Pn+P)/2;


   float t1 = (pe*sex+pn*snx+pw*swx+ps*ssx)*dt/(ro*v);
   U = Ut-(pe*sex+pn*snx+pw*swx+ps*ssx)*dt/(ro*v);
   V = Vt-(pe*sey+pn*sny+pw*swy+ps*ssy)*dt/(ro*v);
}
void vol::AddBToMatr()
   {
   float swy = (a2.x-a1.x);
   float ssy = -(a3.x-a1.x);
   float sny = (a4.x-a2.x);
   float sey = -(a4.x-a3.x);

   float swx = -(a2.y-a1.y);
   float ssx = (a3.y-a1.y);
   float snx = -(a4.y-a2.y);
   float sex = (a4.y-a3.y);
   
      float v = 0.5*((a4.x-a1.x)*(a2.y-a3.y)-(a2.x-a3.x)*(a4.y-a1.y));
      
      float Ue = 0;
      float Un = 0;
      float Uw = 0;
      float Us = 0;
      float Ve = 0;
      float Vn = 0;
      float Vw = 0;
      float Vs = 0;

      Uw = vec[i-1][j].Ut;
      Vw = vec[i-1][j].Vt;
      Us = vec[i][j-1].Ut;
      Vs = vec[i][j-1].Vt;
      Ue = vec[i+1][j].Ut;
      Ve = vec[i+1][j].Vt;
      Un = vec[i][j+1].Ut;
      Vn = vec[i][j+1].Vt;

      float gw = (Ut+Uw)*swx+(Vt+Vw)*swy;
      float ge = (Ut+Ue)*sex+(Vt+Ve)*sey;
      float gn = (Ut+Un)*snx+(Vt+Vn)*sny;
      float gs = (Ut+Us)*ssx+(Vt+Vs)*ssy;
      
      if((bool)(state&128))
      {
         if((bool)(state&1))   
            gw = 0;                                 

//на границах градиент давления равен нулю.
         if((bool)(state&2))   
            ge = 0;
         if((bool)(state&4))
            gn = 0;
         if((bool)(state&8))
            gs = 0;
      }
      if(i == 0)
         gw = UB*swx;
      if(i == nx-1)
         ge = UB*sex;
      b.items[i*(ny)+j] = (ro/dt)*0.5*(ge+gn+gw+gs)/v;
   }
void vol::CompTemp()
   {
      //иксовые и игриковые компоненты площади
   float swy = (a2.x-a1.x);
   float ssy = -(a3.x-a1.x);
   float sny = (a4.x-a2.x);
   float sey = -(a4.x-a3.x);

   float swx = -(a2.y-a1.y);
   float ssx = (a3.y-a1.y);
   float snx = -(a4.y-a2.y);
   float sex = (a4.y-a3.y);

   float v = 0.5*((a4.x-a1.x)*(a2.y-a3.y)-(a2.x-a3.x)*(a4.y-a1.y));

   float Ue = 0;
   float Un = 0;
   float Uw = 0;
   float Us = 0;
   float Ve = 0;
   float Vn = 0;
   float Vw = 0;
   float Vs = 0;
      
   float dxe = 0;      //расстояние до соседней ячейки
   float dxn = 0;
   float dxw = 0;
   float dxs = 0;
   
   float dye = 0;
   float dyn = 0;
   float dyw = 0;
   float dys = 0;
      
   Uw = vec[i-1][j].U;
   Vw = vec[i-1][j].V;
   dxw = c.x-vec[i-1][j].c.x;
   dyw = c.y-vec[i-1][j].c.y;

   Us = vec[i][j-1].U;
   Vs = vec[i][j-1].V;
   dxs = c.x-vec[i][j-1].c.x;
   dys = c.y-vec[i][j-1].c.y;

   Ue = vec[i+1][j].U;
   Ve = vec[i+1][j].V;
   dxe = vec[i+1][j].c.x-c.x;
   dye = vec[i+1][j].c.y-c.y;

   Un = vec[i][j+1].U;
   Vn = vec[i][j+1].V;
   dxn = vec[i][j+1].c.x-c.x;
   dyn = vec[i][j+1].c.y-c.y;
   
   if(state&1)                              
   {         
      dxw = 2*c.x-(a1.x+a2.x);
      dyw = 2*c.y-(a1.y+a2.y);
   }
   if(state&2)
   {            
      dxe = (a4.x+a3.x)-2*c.x;
      dye = (a4.y+a3.y)-2*c.y;
   }
   if(state&4)
   {
      dxn = (a4.x+a2.x)-2*c.x;
      dyn = (a4.y+a2.y)-2*c.y;
   }
   if(state&8)
   {
      dxs = 2*c.x-(a1.x+a3.x);
      dys = 2*c.y-(a1.y+a3.y);
   }

   float mn = dxn*dxn+dyn*dyn;
   float mw = dxw*dxw+dyw*dyw;
   float ms = dxs*dxs+dys*dys;
   float me = dxe*dxe+dye*dye;
      
   float ve = 0.25*((V+Ve)*((V+Ve)*sey+(U+Ue)*sex));
   float vn = 0.25*((V+Vn)*((V+Vn)*sny+(U+Un)*snx));
   float vs = 0.25*((V+Vs)*((V+Vs)*ssy+(U+Us)*ssx));
   float vw = 0.25*((V+Vw)*((V+Vw)*swy+(U+Uw)*swx));
      
   float ue = 0.25*((U+Ue)*((V+Ve)*sey+(U+Ue)*sex));
   float un = 0.25*((U+Un)*((V+Vn)*sny+(U+Un)*snx));
   float us = 0.25*((U+Us)*((V+Vs)*ssy+(U+Us)*ssx));
   float uw = 0.25*((U+Uw)*((V+Vw)*swy+(U+Uw)*swx));

   if((bool)(state&128))
   {
      if((bool)(state&1)&(i>0&i<nx-1))   
      {
         vw = 0;
         uw = 0;      
      }
      if((bool)(state&2))   
      {
         ve = 0;
         ue = 0;      
      }
      if((bool)(state&4))
      {
         vn = 0;
         un = 0;      
      }
      if((bool)(state&8))
      {
         vs = 0;
         us = 0;      
      }
   }

   float Ax;
   float Ay;
   Ax = (ue+un+us+uw)/v;            
   Ay = (ve+vn+vs+vw)/v;

   float De = 0;
   float Dw = 0;
   float Ds = 0;
   float Dn = 0;
   if(sex!=0)
      De = (sex*dxe/me);
   if(swx!=0)
      Dw = (swx*dxw/mw);
   if(ssx!=0)
      Ds = (ssx*dxs/ms);
   if(snx!=0)
      Dn = (snx*dxn/mn);
   
   if(sey!=0)
      De += (sey*dye/me);
   if(swy!=0)
      Dw += (swy*dyw/mw);
   if(ssy!=0)
      Ds += (ssy*dys/ms);
   if(sny!=0)
      Dn += (sny*dyn/mn);

   
   if(state&1)                              
   {         
      Uw = -U;
      Vw = -V;
   }
   if(state&2)
   {            
      Ue = -U;
      Ve = -V;
   }
   if(state&4)
   {
      Un = -U;
      Vn = -V;
   }
   if(state&8)
   {
      Us = -U;
      Vs = -V;
   }


   float Dx = u*((U-Uw)*Dw+(Ue-U)*De+(U-Us)*Ds+(Un-U)*Dn)/v;   
   float Dy = u*((V-Vw)*Dw+(Ve-V)*De+(V-Vs)*Ds+(Vn-V)*Dn)/v;
   Ut = U + (Dx-Ax)*dt;
   Vt = V + (Dy-Ay)*dt;

   float gw = 600*Dw;         //Вот здесь! Если не умножать на 600, то метод расходится, для шага по сетке, меньше 1.
   float gs = 600*Ds;
   float gn = 600*Dn;
   float ge = 600*De;
   if((bool)(state&128))
   {
      if((bool)(state&1)&(i>0&i<nx-1))   
         gw = 0;
      if((bool)(state&2))   
         ge = 0;
      if((bool)(state&4))
         gn = 0;
      if((bool)(state&8))
         gs = 0;
   }

   float a = (-gs);         //это я в таком виде записываю матрицу
   float b = (-gw);
   float c = -(gn-gs-gw+ge);
   float d = (ge);
   float e = (gn);
   MRow temp(a,b,c,d,e);
   Matr[i*(ny)+j] = temp;
}
void cont()      //Функция, которая должна все решать
{
   for(int i = 1;i<nx-1;i++)
      for(int j = 1;j<ny-1;j++)
         vec[i][j].CompTemp();
   for(int i = 1;i<nx-1;i++)
      for(int j = 1;j<ny-1;j++)
         vec[i][j].AddBToMatr();
   BCG();
   for(int i = 1;i<nx-1;i++)
      for(int j = 1;j<ny-1;j++)
         vec[i][j].Correct();
}

[Pphantom]

larts в сообщении #861986 писал(а):

Решение расходится при размере ячейки меньше 1.

Условие Куранта выполнено?

P.S. Код, уж простите, жутковат, лезть разбираться страшно...

[larts]

Понимаю, я бы тоже не полез
Выполнено. Это вообще никак не зависит ни от шага по времени, ни от скорости или вязкости. Только сделаешь сетку помельче - и уже давление как-то криво распределяется.

[Pphantom]

larts в сообщении #861998 писал(а):

Это вообще никак не зависит ни от шага по времени, ни от скорости или вязкости. Только сделаешь сетку помельче - и уже давление как-то криво распределяется.

Так ведь условие и завязано на шаг по координате, шаг по времени и скорость, причем в виде $v \tau/h$. Уж очень похожая картина по всем признакам.

[larts]

Если я правильно понимаю, чтоб условие Куранта выполнялось, при уменьшении шага по сетке я должен уменьшить шаг по времени или скорость. Но я могу делать их сколько угодно маленькими - сходимость от этого не зависит
Или, я все никак не могу сообразить, что Вы хотели сказать..

[Pphantom]

larts в сообщении #862093 писал(а):

Если я правильно понимаю, чтоб условие Куранта выполнялось, при уменьшении шага по сетке я должен уменьшить шаг по времени или скорость.

Да.

larts в сообщении #862093 писал(а):

Но я могу делать их сколько угодно маленькими - сходимость от этого не зависит

Из предыдущего описания было похоже, что зависит...