流体体积法（Volume ofFluid）是一种典型的界面追踪方法，这种方法选取流体体积分数为界面函数S。它通过定义一个体积分数$ C $(指定的流体体积分数占网格体积的百分比)来描述界面。因此只有所在网格体积分数来描述 $ 0<C<1 $ ，而界面两侧的网格内体积分数分别满足 $ C=0 $ 和$ C=1 $。体积分数 $ C $的输运方程为： $$ {\partial C \over \partial t}+ \vec v \cdot \nabla C =0 $$
由体积分数的物理意义可知，可以采取精确的数值算法来构造体积分数的对流量，这样就可以保证VOF方法有很好的守恒性。VOF 方法中对满网格和空网格是容易的，关键问题是够着界面所在的网格即$ 0<C<1 $ 处的数值通量。我们采用两次界面推进两次界面重构的方法，分两次在x 和 y 方向进行推进。
不同网格下计算的结果：
1. 50X50
2. 100 X 100 网格
3. 200 X100 网格
4. 不同网格情况下，一定时间步长时质量损失图。
分析：
可以看出，VOF方法对界面的捕捉能力较好，通过加密网格可以有效的减少质量的损失。
附代码：
其中类和头文件的定义如下：
详细的代码托管在 githup 上
```C++
#include <iostream>
#include <vector>
#include <fstream>
#include <stdlib.h>
#include <string>
#include <sstream> using namespace std;
class vof;
class node;
class element;
extern bool judgev; double check_norline(const double nnx, const double ny, const double alpha, const double pfi, const double dx); #ifndef VOF_H
#define VOF_H
class vof
{
public:
vof();
virtual ~vof(); public:
void initial_nodes();
void calculate(const int & maxTime,int & iouput); private:
void output(int &filenum);
void advance(int & now);
void vof::check_pfi(ofstream & fcout);
void vof::check_pfi(double & fcout); private:
const double PI = 3.14159;
double xmax; //the length of calculate area
double ymax; //the wifth of calculate area
int l; //the node number in x director
int m; //the node number in y director
int ioutput; double dx;
double dy;
double dt;
vector<vector<node> > nodes;
vector<vector<element> > nelem;
}; #endif // VOF_H #ifndef ELEMENT_H
#define ELEMENT_H
class element
{
friend void vof::output(int &);
friend void vof::check_pfi(ofstream & fcout);
friend void vof::check_pfi(double & fcout); public:
element();
virtual ~element();
element& operator=(const element &rhs); public:
void init_pfi(int i, int j, const double dx, vector<vector<node> > & nodes);
void get_norline(const int & ii, const int & jj, const double & dx, vector<vector<element> > & nelem);
double element::get_alpha(const double &nx, const double & ny, const double &pfi, const double &dx);
void ele_xflux(const double &dt, const int &ii, const int &jj, const double &dx, vector<vector<node> > & nodes, int &inow);
void ele_xadvance(const int ii, const int jj, vector<vector<element> > & nelem);
void ele_yflux(const double &dt, const int &ii, const int &jj, const double &dx, vector<vector<node> > & nodes,int &inow);
void ele_yadvance(const int ii, const int jj, vector<vector<element> > & nelem); private:
double element::hea_func(const double & x);
double element::dichotomy_func(const double & aalpha, const double & area, const double & h, const double &nnx, const double& nny);
double element::sign(const double & x);
void element::crossline(const double &x1, const double y1, const double &x2, const double &y2, double &x0, double &y0); private:
const double EPS = 1.0E-9;
const double EPSpfi = 5.0E-3;
double pfi;
double nx;
double ny;
double alpha; double left_flux;
double right_flux; //vector<node> elem_node; //the four nodes of this element }; #endif // ELEMENT_H #ifndef NODE_H
#define NODE_H class node
{
friend void vof::output(int &);
friend void element::ele_xflux(const double &dt, const int &ii, const int &jj, const double &dx, vector<vector<node> > & nodes, int &inow);
friend void element::ele_yflux(const double &dt, const int &ii, const int &jj, const double &dy, vector<vector<node> > & nodes,int &inow);
public:
node(const double &xx, const double & yy);
node& operator=(const node &rhs);
node();
virtual ~node(); public:
double distance();
void initial_node(double x,double y); private:
void get_velicity(); private:
const double PI = 3.1415926;
double x;
double y;
double vx, vy; }; #endif // NODE_H ```