tadd reynolds number function - 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
---
(DIR) commit 0b6294bb9552807ee7c7671c2ca4a073b9f13076
(DIR) parent 78d454149f22e5d2eb8b39a84159b78ff5259dc0
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Fri, 12 Sep 2014 12:38:15 +0200
add reynolds number function
Diffstat:
M python/diffusivity-results.py | 12 +++++++-----
M python/permeability-results.py | 46 +++++++++++++++++++++++++------
M python/shear-results-forces.py | 6 +++---
M python/shear-results.py | 8 +++++---
M python/sphere.py | 18 ++++++++++++++++++
5 files changed, 71 insertions(+), 19 deletions(-)
---
(DIR) diff --git a/python/diffusivity-results.py b/python/diffusivity-results.py
t@@ -53,12 +53,14 @@ ax1.set_ylabel('Hydraulic diffusivity $\\alpha$ [m$^2$s$^{-1}$]')
#ax1.grid()
ax2 = ax1.twinx()
-color = 'b'
-ax2.plot(load, phi_bar, 'o--' + color)
-ax2.set_ylabel('Mean porosity $\\bar{\\phi}$ [-]', color=color)
+#color = 'black'
+#ax2.plot(load, phi_bar, 'o--' + color)
+ax2.plot(load, phi_bar, 'o--', color='black')
+ax2.set_ylabel('Mean porosity $\\bar{\\phi}$ [-]')
+#ax2.set_ylabel('Mean porosity $\\bar{\\phi}$ [-]', color=color)
ax2.get_yaxis().get_major_formatter().set_useOffset(False)
-for tl in ax2.get_yticklabels():
- tl.set_color(color)
+#for tl in ax2.get_yticklabels():
+ #tl.set_color(color)
filename = 'diffusivity-sigma0-vs-alpha.pdf'
plt.tight_layout()
(DIR) diff --git a/python/permeability-results.py b/python/permeability-results.py
t@@ -2,6 +2,8 @@
import matplotlib
matplotlib.use('Agg')
matplotlib.rcParams.update({'font.size': 18, 'font.family': 'serif'})
+matplotlib.rc('text', usetex=True)
+matplotlib.rcParams['text.latex.preamble']=[r"\usepackage{amsmath}"]
import shutil
import os
t@@ -18,6 +20,7 @@ K = [[], [], []]
dpdz = [[], [], []]
Q = [[], [], []]
phi_bar = [[], [], []]
+Re = [[], [], []]
c = 0
t@@ -35,6 +38,7 @@ for c_grad_p in cvals:
dpdz[c] = numpy.zeros_like(K[c])
Q[c] = numpy.zeros_like(K[c])
phi_bar[c] = numpy.zeros_like(K[c])
+ Re[c] = numpy.zeros_like(K[c])
i = 0
for sid in sids:
t@@ -49,6 +53,10 @@ for c_grad_p in cvals:
pc.findMeanPorosity()
#pc.plotEvolution()
phi_bar[c][i] = pc.phi_bar
+
+ sim = sphere.sim(sid, fluid=True)
+ sim.readlast(verbose=False)
+ Re[c][i] = numpy.mean(sim.ReynoldsNumber())
else:
print(sid + ' not found')
t@@ -60,19 +68,38 @@ for c_grad_p in cvals:
#sim.writeVTKall()
c += 1
-fig = plt.figure()
+fig = plt.figure(figsize=(8,12))
-#plt.subplot(3,1,1)
-plt.xlabel('Pressure gradient $\\Delta p/\\Delta z$ [kPa m$^{-1}$]')
-plt.ylabel('Hydraulic conductivity $K$ [ms$^{-1}$]')
-plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+ax1 = plt.subplot(3,1,1)
+ax2 = plt.subplot(3,1,2, sharex=ax1)
+ax3 = plt.subplot(3,1,3, sharex=ax1)
+lines = ['-', '--', '-.', ':']
+markers = ['o', 'x', '^', '+']
for c in range(len(cvals)):
dpdz[c] /= 1000.0
#plt.plot(dpdz[c], K[c], 'o-', label='$c$ = %.2f' % (cvals[c]))
#plt.semilogx(dpdz[c], K[c], 'o-', label='$c$ = %.2f' % (cvals[c]))
#plt.semilogy(dpdz[c], K[c], 'o-', label='$c$ = %.2f' % (cvals[c]))
- plt.loglog(dpdz[c], K[c], 'o-', label='$c$ = %.2f' % (cvals[c]))
-plt.grid()
+ ax1.loglog(dpdz[c], K[c], label='$c$ = %.2f' % (cvals[c]),
+ linestyle=lines[c], marker=markers[c], color='black')
+ ax2.semilogx(dpdz[c], phi_bar[c], label='$c$ = %.2f' % (cvals[c]),
+ linestyle=lines[c], marker=markers[c], color='black')
+ ax3.loglog(dpdz[c], Re[c], label='$c$ = %.2f' % (cvals[c]),
+ linestyle=lines[c], marker=markers[c], color='black')
+
+ax1.set_ylabel('Hydraulic conductivity $K$ [ms$^{-1}$]')
+#ax1.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+
+ax2.set_ylabel('Mean porosity $\\bar{\\phi}$ [-]')
+
+ax3.set_ylabel('Mean Reynolds number $\\bar{Re}$ [-]')
+
+ax3.set_xlabel('Pressure gradient $\\Delta p/\\Delta z$ [kPa m$^{-1}$]')
+plt.setp(ax1.get_xticklabels(), visible=False)
+
+ax1.grid()
+ax2.grid()
+ax3.grid()
#plt.subplot(3,1,2)
#plt.xlabel('Pressure gradient $\\Delta p/\\Delta z$ [Pa m$^{-1}$]')
t@@ -86,7 +113,10 @@ plt.grid()
#plt.plot(dpdz, phi_bar, '+')
#plt.grid()
-plt.legend(loc='lower left', prop={'size':18}, fancybox=True, framealpha=0.5)
+ax1.legend(loc='best', prop={'size':18}, fancybox=True, framealpha=0.5)
+ax2.legend(loc='best', prop={'size':18}, fancybox=True, framealpha=0.5)
+ax3.legend(loc='best', prop={'size':18}, fancybox=True, framealpha=0.5)
+
plt.tight_layout()
filename = 'permeability-dpdz-vs-K-vs-c.pdf'
#print(os.getcwd() + '/' + filename)
(DIR) diff --git a/python/shear-results-forces.py b/python/shear-results-forces.py
t@@ -15,10 +15,10 @@ from matplotlib.ticker import MaxNLocator
#steps = [5, 10, 100]
#steps = [5, 10]
-steps = sys.argv[1:]
+steps = sys.argv[2:]
nsteps_avg = 5 # no. of steps to average over
-sigma0 = 10.0e3
+sigma0 = float(sys.argv[1])
c_grad_p = 1.0
c_phi = 1.0
t@@ -179,7 +179,7 @@ plt.tight_layout()
plt.subplots_adjust(wspace = .05)
plt.MaxNLocator(nbins=4)
-filename = 'shear-10kPa-forces.pdf'
+filename = 'shear-' + str(int(sigma0/1000.0)) + 'kPa-forces.pdf'
plt.savefig(filename)
shutil.copyfile(filename, '/home/adc/articles/own/2-org/' + filename)
print(filename)
(DIR) diff --git a/python/shear-results.py b/python/shear-results.py
t@@ -7,13 +7,15 @@ matplotlib.rcParams['text.latex.preamble']=[r"\usepackage{amsmath}"]
import shutil
import os
+import sys
import numpy
import sphere
from permeabilitycalculator import *
import matplotlib.pyplot as plt
#sigma0_list = numpy.array([1.0e3, 2.0e3, 4.0e3, 10.0e3, 20.0e3, 40.0e3])
-sigma0 = 10.0e3
+#sigma0 = 10.0e3
+sigma0 = float(sys.argv[1])
#cvals = [1.0, 0.1]
cvals = [1.0]
c_phi = 1.0
t@@ -28,7 +30,7 @@ p_max = [[], [], []]
fluid=True
# dry shear
-sid = 'shear-sigma0=' + str(10.0e3)
+sid = 'shear-sigma0=' + sys.argv[1]
sim = sphere.sim(sid)
sim.readlast(verbose=False)
sim.visualize('shear')
t@@ -129,7 +131,7 @@ ax3.legend(loc='lower right', prop={'size':18}, fancybox=True,
framealpha=legend_alpha)
plt.tight_layout()
-filename = 'shear-10kPa-stress-dilation.pdf'
+filename = 'shear-' + str(int(sigma0/1000.0)) + 'kPa-stress-dilation.pdf'
#print(os.getcwd() + '/' + filename)
plt.savefig(filename)
shutil.copyfile(filename, '/home/adc/articles/own/2-org/' + filename)
(DIR) diff --git a/python/sphere.py b/python/sphere.py
t@@ -4387,6 +4387,24 @@ class sim:
plt.clf()
plt.close(fig)
+ def ReynoldsNumber(self):
+ '''
+ Estimate the per-cell Reynolds number by: Re = rho * ||v_f|| * dx/mu
+
+ :returns: Reynolds number
+ :return type: Numpy array with dimensions like the fluid grid
+ '''
+
+ # find magnitude of fluid velocity vectors
+ self.v_f_magn = numpy.empty_like(self.p_f)
+ for z in numpy.arange(self.num[2]):
+ for y in numpy.arange(self.num[1]):
+ for x in numpy.arange(self.num[0]):
+ self.v_f_magn = self.v_f[x,y,z,:].dot(self.v_f[x,y,z,:])
+
+ self.Re = self.rho_f*self.v_f_magn*self.L[0]/self.num[0]/self.mu
+ return self.Re
+
def plotLoadCurve(self, graphics_format='png'):
'''
Plot the load curve (log time vs. upper wall movement). The plot is