VOF 方法捕捉界面--粘性剪切流动算例

流体体积法（Volume ofFluid）是一种典型的界面追踪方法，这种方法选取流体体积分数为界面函数S。它通过定义一个体积分数$ C $(指定的流体体积分数占网格体积的百分比)来描述界面。因此只有所在网格体积分数来描述 $ 0<C<1 $ ，而界面两侧的网格内体积分数分别满足 $ C=0 $ 和$ C=1 $。体积分数 $ C $的输运方程为： $$ {partial C over partial t}+ vec v cdot abla C =0 $$

由体积分数的物理意义可知，可以采取精确的数值算法来构造体积分数的对流量，这样就可以保证VOF方法有很好的守恒性。VOF 方法中对满网格和空网格是容易的，关键问题是够着界面所在的网格即$ 0<C<1 $ 处的数值通量。我们采用两次界面推进两次界面重构的方法，分两次在x 和 y 方向进行推进。

不同网格下计算的结果：

1. 50X50

2. 100 X 100 网格

3. 200 X100 网格

4. 不同网格情况下，一定时间步长时质量损失图。

分析：

可以看出，VOF方法对界面的捕捉能力较好，通过加密网格可以有效的减少质量的损失。

附代码：

其中类和头文件的定义如下：

详细的代码托管在 githup 上

  1 ```C++
  2 #include <iostream>
  3 #include <vector>
  4 #include <fstream>
  5 #include <stdlib.h>
  6 #include <string>
  7 #include <sstream>
  8 
  9 using namespace std;
 10 class vof;
 11 class node;
 12 class element;
 13 extern  bool judgev;
 14 
 15 double  check_norline(const double nnx, const double ny, const double  alpha, const double  pfi, const double dx);
 16 
 17 #ifndef VOF_H
 18 #define VOF_H
 19 class vof
 20 {
 21 public:
 22     vof();
 23     virtual ~vof();
 24 
 25 public:
 26     void initial_nodes();
 27     void calculate(const int  & maxTime,int & iouput);
 28 
 29 private:
 30     void output(int &filenum);
 31     void advance(int & now);
 32     void vof::check_pfi(ofstream & fcout);
 33     void vof::check_pfi(double & fcout);
 34 
 35 private:     
 36     const double PI = 3.14159;
 37     double xmax;    //the length of  calculate area
 38     double ymax;    //the wifth of calculate area
 39     int l;      //the node number in x director
 40     int m;      //the node number in y director
 41     int ioutput;
 42 
 43     double dx;
 44     double dy;
 45     double dt;
 46     vector<vector<node> > nodes;
 47     vector<vector<element> > nelem;
 48 };
 49 
 50 #endif // VOF_H
 51 
 52 #ifndef ELEMENT_H
 53 #define ELEMENT_H
 54 class element
 55 {
 56     friend void vof::output(int &);
 57     friend void vof::check_pfi(ofstream & fcout);
 58     friend void vof::check_pfi(double & fcout);
 59 
 60 public:
 61     element();
 62     virtual ~element();
 63     element& operator=(const element &rhs);
 64 
 65 public:
 66     void init_pfi(int i, int j, const double dx, vector<vector<node> >  & nodes);
 67     void get_norline(const int & ii, const int &  jj, const double & dx, vector<vector<element> > & nelem);
 68     double element::get_alpha(const double &nx, const double & ny, const double &pfi, const double &dx);
 69     void ele_xflux(const double &dt, const int &ii, const int &jj, const double &dx, vector<vector<node> > & nodes, int &inow);
 70     void ele_xadvance(const int ii, const int jj, vector<vector<element> >  & nelem);
 71     void ele_yflux(const double &dt, const int &ii, const int &jj, const double &dx, vector<vector<node> > & nodes,int &inow);
 72     void ele_yadvance(const int ii, const int jj, vector<vector<element> >  & nelem);
 73 
 74 private:
 75     double element::hea_func(const double & x);
 76     double element::dichotomy_func(const double & aalpha, const double & area, const double & h, const double &nnx, const double& nny);
 77     double element::sign(const double & x);
 78     void element::crossline(const double &x1, const double y1, const double &x2, const double &y2, double &x0, double &y0);
 79 
 80 private:
 81     const double EPS = 1.0E-9;
 82     const double EPSpfi = 5.0E-3;
 83     double pfi;
 84     double nx;
 85     double ny;
 86     double alpha;
 87 
 88     double left_flux;
 89     double right_flux;
 90 
 91     //vector<node> elem_node;   //the four nodes of this element
 92 
 93 };
 94 
 95 #endif // ELEMENT_H
 96 
 97 #ifndef NODE_H
 98 #define NODE_H
 99 
100 class node
101 {
102     friend void vof::output(int &);
103     friend void element::ele_xflux(const double &dt, const int &ii, const int &jj, const double &dx, vector<vector<node> > & nodes, int &inow);
104     friend void element::ele_yflux(const double &dt, const int &ii, const int &jj, const double &dy, vector<vector<node> > & nodes,int &inow);
105 public:
106     node(const double &xx, const double & yy);
107     node& operator=(const node &rhs);
108     node();
109     virtual ~node();
110 
111 public:
112     double distance();
113     void initial_node(double x,double y);
114     
115 private:
116     void get_velicity();
117 
118 private:
119     const double PI = 3.1415926;
120     double x;
121     double y;
122     double vx, vy;
123 
124 };
125 
126 #endif // NODE_H
127 
128 ```