流體體積法(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 上
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 ```