tpermeabilitycalculator.py - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
(HTM) git clone git://src.adamsgaard.dk/sphere
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tpermeabilitycalculator.py (6201B)
---
1 #!/usr/bin/env python
2 import sys
3 import sphere
4 import numpy
5 import matplotlib.pyplot as plt
6
7 class PermeabilityCalc:
8 ''' Darcy's law: Q = -k*A/mu * dP '''
9
10 def __init__(self, sid, plot_evolution=True, print_results=True,
11 verbose=True):
12 self.sid = sid
13 self.readfile(verbose)
14 self.findPermeability()
15 self.findConductivity()
16 self.findMeanPorosity()
17 if print_results:
18 self.printResults()
19 if plot_evolution:
20 self.plotEvolution()
21
22 def readfile(self, verbose=True):
23 self.sim = sphere.sim(self.sid, fluid=True)
24 self.sim.readlast(verbose=verbose)
25
26 def findPermeability(self):
27 self.findCellSpacing()
28 self.findCrossSectionalArea()
29 self.findCrossSectionalFlux()
30 self.findPressureGradient()
31 self.k = -self.Q*self.sim.mu/(self.A*self.dPdL) # m^2
32
33 def findConductivity(self):
34 # hydraulic conductivity
35 self.findCellSpacing()
36 self.findCrossSectionalArea()
37 self.findCrossSectionalFlux()
38 self.findPressureGradient()
39 #self.K = self.k/self.sim.mu # m/s
40 self.K = -self.Q * self.dL / (self.A * self.dP)
41
42 def conductivity(self):
43 return self.K[2]
44
45 def c_grad_p(self):
46 return self.sim.c_grad_p[0]
47
48 def findMeanPorosity(self):
49 ''' calculate mean porosity in cells beneath the top wall '''
50
51 if (self.sim.nw > 0):
52 wall_iz = int(self.sim.w_x[0]/self.dx[2])
53 self.phi_bar = numpy.mean(self.sim.phi[:,:,0:wall_iz-1])
54 else:
55 self.phi_bar = numpy.mean(self.sim.phi[:,:,0:-3])
56
57
58 def findCrossSectionalArea(self):
59 ''' Cross sectional area normal to each axis '''
60 self.A = numpy.array([
61 self.sim.L[1]*self.sim.L[2],
62 self.sim.L[0]*self.sim.L[2],
63 self.sim.L[0]*self.sim.L[1]])
64
65 def findCellSpacing(self):
66 self.dx = numpy.array([
67 self.sim.L[0]/self.sim.num[0],
68 self.sim.L[1]/self.sim.num[1],
69 self.sim.L[2]/self.sim.num[2]])
70
71 def findCrossSectionalFlux(self):
72 ''' Flux along each axis, measured at the outer boundaries '''
73 #self.Q = numpy.array([
74 #numpy.mean(self.sim.v_f[-1,:,:]),
75 #numpy.mean(self.sim.v_f[:,-1,:]),
76 #numpy.mean(self.sim.v_f[:,:,-1])])*self.A
77
78 self.Q = numpy.zeros(3)
79
80 self.A_cell = numpy.array([
81 self.dx[1]*self.dx[2],
82 self.dx[0]*self.dx[2],
83 self.dx[0]*self.dx[1]])
84
85 # x axis (0)
86 for y in numpy.arange(self.sim.num[1]):
87 for z in numpy.arange(self.sim.num[2]):
88 self.Q[0] += self.sim.v_f[-1,y,z,0] * self.A_cell[0]
89
90 # y axis (1)
91 for x in numpy.arange(self.sim.num[0]):
92 for z in numpy.arange(self.sim.num[2]):
93 self.Q[1] += self.sim.v_f[x,-1,z,1] * self.A_cell[1]
94
95 # z axis (2)
96 for x in numpy.arange(self.sim.num[0]):
97 for y in numpy.arange(self.sim.num[1]):
98 self.Q[2] += self.sim.v_f[x,y,-1,2] * self.A_cell[2]
99
100 def findPressureGradient(self):
101 ''' Determine pressure gradient by finite differencing the
102 mean values at the outer boundaries '''
103 self.dP = numpy.array([
104 numpy.mean(self.sim.p_f[-1,:,:]) - numpy.mean(self.sim.p_f[0,:,:]),
105 numpy.mean(self.sim.p_f[:,-1,:]) - numpy.mean(self.sim.p_f[:,0,:]),
106 numpy.mean(self.sim.p_f[:,:,-1]) - numpy.mean(self.sim.p_f[:,:,0])
107 ])
108 self.dL = self.sim.L
109 self.dPdL = self.dP/self.dL
110
111 def printResults(self):
112 print('\n### Permeability resuts for "' + self.sid + '" ###')
113 print('Pressure gradient: dPdL = ' + str(self.dPdL) + ' Pa/m')
114 print('Flux: Q = ' + str(self.Q) + ' m^3/s')
115 print('Intrinsic permeability: k = ' + str(self.k) + ' m^2')
116 print('Saturated hydraulic conductivity: K = ' + str(self.K) + ' m/s')
117 print('Mean porosity: phi_bar = ' + str(self.phi_bar) + '\n')
118
119 def plotEvolution(self, axis=2, outformat='png'):
120 '''
121 Plot temporal evolution of parameters on the selected axis.
122 Note that the first 5 output files are ignored.
123 '''
124 skipsteps = 5
125 nsteps = self.sim.status() - skipsteps
126 self.t_series = numpy.empty(nsteps)
127 self.Q_series = numpy.empty((nsteps, 3))
128 self.phi_bar_series = numpy.empty(nsteps)
129 self.k_series = numpy.empty((nsteps, 3))
130 self.K_series = numpy.empty((nsteps, 3))
131
132 print('Reading ' + str(nsteps) + ' output files... '),
133 sys.stdout.flush()
134 for i in numpy.arange(skipsteps, self.sim.status()):
135 self.sim.readstep(i, verbose=False)
136
137 self.t_series[i-skipsteps] = self.sim.time_current[0]
138
139 self.findCrossSectionalFlux()
140 self.Q_series[i-skipsteps,:] = self.Q
141
142 self.findMeanPorosity()
143 self.phi_bar_series[i-skipsteps] = self.phi_bar
144
145 self.findPermeability()
146 self.k_series[i-skipsteps,:] = self.k
147
148 self.findConductivity()
149 self.K_series[i-skipsteps,:] = self.K
150 print('Done')
151
152 fig = plt.figure()
153
154 plt.subplot(2,2,1)
155 plt.xlabel('Time $t$ [s]')
156 plt.ylabel('Flux $Q$ [m^3/s]')
157 plt.plot(self.t_series, self.Q_series[:,axis])
158 #plt.legend()
159 plt.grid()
160
161 plt.subplot(2,2,2)
162 plt.xlabel('Time $t$ [s]')
163 plt.ylabel('Porosity $\phi$ [-]')
164 plt.plot(self.t_series, self.phi_bar_series)
165 plt.grid()
166
167 plt.subplot(2,2,3)
168 plt.xlabel('Time $t$ [s]')
169 plt.ylabel('Permeability $k$ [m^2]')
170 plt.plot(self.t_series, self.k_series[:,axis])
171 plt.grid()
172
173 plt.subplot(2,2,4)
174 plt.xlabel('Time $t$ [s]')
175 plt.ylabel('Conductivity $K$ [m/s]')
176 plt.plot(self.t_series, self.K_series[:,axis])
177 plt.grid()
178
179 fig.tight_layout()
180
181 filename = self.sid + '-permeability.' + outformat
182 plt.savefig(filename)
183 print('Figure saved as "' + filename + '"')