Main author: Wim Degruyter
This tutorial illustrates, step-by-step, how to compute the permeability of a given porous media. This process consists of three major steps:
The code developed in this tutorial can be found in the directory palabos/examples/tutorial/permeability. This Palabos code has been developed as part of a permeability study of vesicular volcanic rocks [1].
[1] | Degruyter W., Bachmann O., Burgisser A., Malaspinas O. A synchrotron perspective on gas flow through pumices. submitted to Geosphere |
---|
We assume you already have downloaded and installed Palabos. Installation support is provided in the user’s guide. We also recommend to go through the other Palabos tutorials to get acquainted with compiling and running of Palabos programs.
The input required for permeability are a series of black and white images in .bmp format, with filename namexxxx.bmp with xxxx numbers starting from 0001. In the directory twoSpheres one can find a series of bw images defining two connected hemispheres. To inspect the stack of images one can use e.g. ImageJ. Here are two example slices, inlet (left) and halfway (right):
Convert image stack to .dat file
The first step is to create an input file from the images, so the code is able to read in the geometry. This conversion is done by a Matlab script called createDAT.m. Open Matlab and go to the directory of createDAT.m. At the prompt type createDAT(number-of-files, path/to/inputprefix,path/to/output.dat) e.g. try for the twoSpheres files and type createDAT(48, 'twoSpheres/', 'twoSpheres', 'twoSpheres.dat'). Every voxel is given a value: 0 for a fluid voxel (blue), 1 for a material voxel that touches (26-connected) a pore voxel (green), and 2 for an interior material voxel (red) illustrated again by the inlet slice (left) and the halfway slice (right):
The Matlab script visualizes the conversion fo illustration purposes. If you are converting a large number of files this process can take some time and it is better to switch of the display by commenting out lines 113-115, 258-260, and 372-374.
Once the .dat file is created we can start a simulation.
The permeability.cpp code is listed below:
1. #include "palabos3D.h"
2. #include "palabos3D.hh"
3. #include <vector>
4. #include <cmath>
5.
6. using namespace plb;
7. using namespace std;
8.
9.
10. typedef double T;
11. #define DESCRIPTOR descriptors::D3Q19Descriptor
12.
13.
14. // This function object returns a zero velocity, and a pressure which decreases
15. // linearly in x-direction. It is used to initialize the particle populations.
16. class PressureGradient {
17. public:
18. PressureGradient(T deltaP_, plint nx_) : deltaP(deltaP_), nx(nx_)
19. { }
20. void operator() (plint iX, plint iY, plint iZ, T& density, Array<T,3>& velocity) const
21. {
22. velocity.resetToZero();
23. density = 1. - deltaP*DESCRIPTOR<T>::invCs2 / (T)(nx-1) * (T)iX;
24.
25. }
26. private:
27. T deltaP;
28. plint nx;
29. };
30.
31. void porousMediaSetup( MultiBlockLattice3D<T,DESCRIPTOR>& lattice,
32. OnLatticeBoundaryCondition3D<T,DESCRIPTOR>* boundaryCondition,
33. MultiScalarField3D<int>& geometry, T deltaP)
34. {
35. const plint nx = lattice.getNx();
36. const plint ny = lattice.getNy();
37. const plint nz = lattice.getNz();
38.
39. pcout << "Definition of inlet/outlet." << endl;
40. Box3D inlet (0,0, 1,ny-2, 1,nz-2);
41. boundaryCondition->addPressureBoundary0N(inlet, lattice);
42. setBoundaryDensity(lattice, inlet, 1.);
43.
44. Box3D outlet(nx-1,nx-1, 1,ny-2, 1,nz-2);
45. boundaryCondition->addPressureBoundary0P(outlet, lattice);
46. setBoundaryDensity(lattice, outlet, 1. - deltaP*DESCRIPTOR<T>::invCs2);
47.
48. pcout << "Definition of the geometry." << endl;
49. // Where "geometry" evaluates to 1, use bounce-back.
50. defineDynamics(lattice, geometry, new BounceBack<T,DESCRIPTOR>(), 1);
51. // Where "geometry" evaluates to 2, use no-dynamics (which does nothing).
52. defineDynamics(lattice, geometry, new NoDynamics<T,DESCRIPTOR>(), 2);
53.
54. pcout << "Initilization of rho and u." << endl;
55. initializeAtEquilibrium( lattice, lattice.getBoundingBox(), PressureGradient(deltaP, nx) );
56.
57. lattice.initialize();
58. delete boundaryCondition;
59. }
60.
61. void writeGifs(MultiBlockLattice3D<T,DESCRIPTOR>& lattice, plint iter)
62. {
63. const plint nx = lattice.getNx();
64. const plint ny = lattice.getNy();
65. const plint nz = lattice.getNz();
66.
67. const plint imSize = 600;
68. ImageWriter<T> imageWriter("leeloo");
69.
70. // Write velocity-norm at x=0.
71. imageWriter.writeScaledGif( createFileName("ux_inlet", iter, 6),
72. *computeVelocityNorm(lattice, Box3D(0,0, 0,ny-1, 0,nz-1)),
73. imSize, imSize );
74.
75. // Write velocity-norm at x=nx/2.
76. imageWriter.writeScaledGif( createFileName("ux_half", iter, 6),
77. *computeVelocityNorm(lattice, Box3D(nx/2,nx/2, 0,ny-1, 0,nz-1)),
78. imSize, imSize );
79.
80. }
81.
82.
83. void writeVTK(MultiBlockLattice3D<T,DESCRIPTOR>& lattice, plint iter)
84. {
85. VtkImageOutput3D<T> vtkOut(createFileName("vtk", iter, 6), 1.);
86. vtkOut.writeData<float>(*computeVelocityNorm(lattice), "velocityNorm", 1.);
87. vtkOut.writeData<3,float>(*computeVelocity(lattice), "velocity", 1.);
88. }
89.
90. T computePermeability (
91. MultiBlockLattice3D<T,DESCRIPTOR>& lattice, T nu, T deltaP, Box3D domain )
92. {
93. pcout << "Computing the permeability." << endl;
94.
95. // Compute only the x-direction of the velocity (direction of the flow).
96. plint xComponent = 0;
97. plint nx = lattice.getNx();
98.
99. T meanU = computeAverage (
100. *computeVelocityComponent (lattice, domain, xComponent ) );
101.
102.
103. pcout << "Average velocity = " << meanU << endl;
104. pcout << "Lattice viscosity nu = " << nu << endl;
105. pcout << "Grad P = " << deltaP/(T)(nx-1) << endl;
106. pcout << "Permeability = " << nu*meanU / (deltaP/(T)(nx-1)) << endl;
107.
108. return meanU;
109. }
110.
111. int main(int argc, char **argv)
112. {
113. plbInit(&argc, &argv);
114.
115. if (argc!=7)
116. {
117. pcout << "Error missing some input parameter\n";
118. pcout << "The structure is :\n";
119. pcout << "1. Input file name.\n";
120. pcout << "2. Output directory name.\n";
121. pcout << "3. number of cells in X direction.\n";
122. pcout << "4. number of cells in Y direction.\n";
123. pcout << "5. number of cells in Z direction.\n";
124. pcout << "6. Delta P .\n";
125. exit (EXIT_FAILURE);
126. }
127. std::string fNameIn = argv[1];
128. std::string fNameOut = argv[2];
129.
130. const plint nx = atoi(argv[3]);
131. const plint ny = atoi(argv[4]);
132. const plint nz = atoi(argv[5]);
133. const T deltaP = atof(argv[6]);
134.
135. global::directories().setOutputDir(fNameOut+"/");
136.
137. const T omega = 1.0;
138. const T nu = ((T)1/omega-0.5)/DESCRIPTOR<T>::invCs2;
139.
140. pcout << "Creation of the lattice." << endl;
141. MultiBlockLattice3D<T,DESCRIPTOR> lattice(nx,ny,nz, new BGKdynamics<T,DESCRIPTOR>(omega));
142. // Switch off periodicity.
143. lattice.periodicity().toggleAll(false);
144.
145. pcout << "Reading the geometry file." << endl;
146. MultiScalarField3D<int> geometry(nx,ny,nz);
147. plb_ifstream geometryFile(fNameIn.c_str());
148. if (!geometryFile.is_open()) {
149. pcout << "Error: could not open geometry file " << fNameIn << endl;
150. return -1;
151. }
152. geometryFile >> geometry;
153.
154. pcout << "nu = " << nu << endl;
155. pcout << "deltaP = " << deltaP << endl;
156. pcout << "omega = " << omega << endl;
157. pcout << "nx = " << lattice.getNx() << endl;
158. pcout << "ny = " << lattice.getNy() << endl;
159. pcout << "nz = " << lattice.getNz() << endl;
160.
161. porousMediaSetup(lattice, createLocalBoundaryCondition3D<T,DESCRIPTOR>(), geometry, deltaP);
162.
163. // The value-tracer is used to stop the simulation once it has converged.
164. // 1st parameter:velocity
165. // 2nd parameter:size
166. // 3rd parameters:threshold
167. // 1st and second parameters ae used for the length of the time average (size/velocity)
168. util::ValueTracer<T> converge(1.0,1000.0,1.0e-4);
169.
170. pcout << "Simulation begins" << endl;
171. plint iT=0;
172.
173. const plint maxT = 30000;
174. for (;iT<maxT; ++iT)
175. {
176. if (iT % 20 == 0) {
177. pcout << "Iteration " << iT << endl;
178. }
179. if (iT % 500 == 0 && iT>0) {
180. writeGifs(lattice,iT);
181. }
182.
183. lattice.collideAndStream();
184. converge.takeValue(getStoredAverageEnergy(lattice),true);
185.
186. if (converge.hasConverged())
187. {
188. break;
189. }
190. }
191.
192. pcout << "End of simulation at iteration " << iT << endl;
193.
194. pcout << "Permeability:" << endl << endl;
195. computePermeability(lattice, nu, deltaP, lattice.getBoundingBox());
196. pcout << endl;
197.
198. pcout << "Writing VTK file ..." << endl << endl;
199. writeVTK(lattice,iT);
200. pcout << "Finished!" << endl << endl;
201. }
Line 1-12General definitions and includes are discussed in the other tutorials of Palabos. Here we use the D3Q19 lattice, but other lattices (e.g. D3Q15 or D3Q27) are also possible.Line 16-29This functional defines the initial conditions. It is used in function porousMediaSetup: the fluid has initially zero velocity, with a linear pressure gradient in the x-direction.Line 31-59Here, the boundary conditions are set, i.e. all voxels read in as 1 are defined to have bounce back boundary conditions, and all voxels read in with value two are set to carry no dynamics. The inlet and outlet are defined to have a fixed pressure difference.Line 61-109Various outputs are created. During the simulation .gif files showing the velocity distribution within a slice are made (line 61-80). A .vti file with the velocity distribution is written, which can be read by Paraview (line 83-88). The permeability is computed by applying Darcy’s law to the simulated velocity data (line 90-109).Line 111-201Main part of the code. The code requires 6 input values from the user: the input.dat file, the output directory, the size of the file in each orthogonal drection, and finally the pressure difference between the in- and outlet. From these input values, the whole run is instantiated and the simulation starts (line 174). The ValueTracer (line 168) is used to stop the simulation when steady state is reached. A maximum of number of iterations is defined (line 173). Every 500 steps a .gif file is written to the output directory (line 179). Once converged the permeability of the medium is written to the screen and a .vti file is created.
First the code needs to be compiled. Check if the Makefile is in order. Open a terminal and type make, once in the permeability directory. To run a simulation type ./permeability path/to/input.dat path/to/outputdir/ nx ny nz deltaP in a terminal in the permeability directory. Let’s test it on our example by typing ./permeability twoSpheres.dat tmp/ 48 64 64 0.00005. The progress of the simulation is written to the screen. One can find the .gif and .vti files in the tmp subdirectory.
Use our favorite image program to visualize the .gif files. Here we see the velocity distribution at the inlet (left) and halfway (right) after 4500 iterations:
Paraview allows you to visualize the .vti file in 3D. Here are some examples:
本文分享自 传输过程数值模拟学习笔记 微信公众号,前往查看
如有侵权,请联系 cloudcommunity@tencent.com 删除。
本文参与 腾讯云自媒体同步曝光计划 ,欢迎热爱写作的你一起参与!