流体体积法（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 ```