import pandas as pd
from os.path import join
import numpy as np
from os import listdir
path = 'data/grainsize_distributions'
sites = listdir(path)
sites = [s for s in sites if not s.endswith('.csv')]
column_names = ['d_lower','vol_frac','-2SD','+2SD','d_center',\
'd_upper','d_center_phi','cum']
datasets = pd.DataFrame()
for site in sites:
site_path = join(path, site)
data_files = listdir(site_path)
data_files.sort()
for df in data_files:
# z-position of the sample for which the grainsize
# distribution was measured (this is an integer number,
# NOT a depth in centimeters)
z = int(df.split(',')[1].replace(').csv', ''))
data = pd.read_csv(join(site_path, df), header=[0,1,2,3])
data.columns = column_names
data['z'] = z
data['site'] = site
# convert data to micrometers
data['d_center'] = data['d_center']*1e-6
data['d_upper'] = data['d_upper']*1e-6
data['d_lower'] = data['d_lower']*1e-6
# calculate area fraction from volume fraction
data['radius'] = data['d_center'] / 2
area_frac = 3 * data['vol_frac'] / data['radius']
area_frac = area_frac / area_frac.sum()
data['area_frac'] = area_frac
datasets = datasets.append(data, ignore_index=True)
From [1], p.47: Given a non-spherical particle with some volume V and some surface area $A$, $d_S$ is defined as the diameter of the sphere with the same surface to area ratio than the original particle. This concept can be transferred to a porous medium with a distribution of particle diameters: Given a polydisperse porous medium with total particle volume $V$ total and particle surface area A total and some distribution of particle diameters $d_i$. The Sauter diameter of the medium is defined as the particle diameter $d_S$ of a monodisperse porous medium that gives the same ratio of total volume to total surface area: \begin{align} d_S := 6\cdot \frac{V_{total}}{A_{total}}\,. \end{align} To calculate the $d_S$ from LPS data, we first need to calculate the total volume of the sample. We do this by calculating the number $n_i$ of particles in each bin $i$ \begin{align} n_i = \varphi_i \cdot\frac{V_{total}}{V_i}\,. \end{align} This is a system of $N$ equations which is closed because $\sum_{i=1}^N\varphi_i = 1$ by definition. We solve for $V_{total}$ by arbitrarily assuming that $n_N = 1$ in the highest bin that recorded data. This is possible because we are looking for a ratio rather than an absolute value. We then calculate the $n_i$ for all bins and subsequently the total area given as \begin{align*} A_{total} = \pi\sum_{i=1}^N n_id_i^2\;. \end{align*} Using the definition of the Sauter diameter above, we can then calculate $d_S$, which gives me a sensible estimate of the intrinsic length scale of the porous medium under the condition that drag is constant.
The main source of error for the calculation of the Sauter diameter comes from the fact that grain sizes are binned by the LPS into bins of varying size and the distribution of grain sizes within one bin is unknown. We estimate error bounds for the Sauter diameter by calculating a sample's Sauter diameter once by assuming all grain diameters are grouped on the low end of each bin and once by assuming the opposite.
def CalculateSauterDiameter(sample_data, diameter_col):
try:
last_chan = sample_data[sample_data['vol_frac'] != 0].index[-1]
except: # if data = nan
return np.nan
d = sample_data[diameter_col] # channel diameter used as grain diameter
d_max = sample_data[diameter_col].loc[last_chan] # last (largest) channel in which grains were recorded
phi_V = sample_data['vol_frac'] # volume fractions in each channel
phi_V_max = sample_data['vol_frac'].loc[last_chan] # volume fraction of the largest recorded grains
# calculate Sauter diameter
const = np.pi / 6.
V = const * d**3 # volume of a sphere with diameter d
V_tot_max = const * d_max**3 / phi_V_max # total volume in the last channel
N = (phi_V * V_tot_max) / (d**3 * const)*1e-2 # number of particles in each bin
V_tot = N * V # total volume
A_tot = N * np.pi * d**2 # total area
d_s = 6 * V_tot.sum() / A_tot.sum() # definition of d_s
return (d_s)
sauter_diameter = pd.DataFrame()
# iterate over all
for site in datasets['site'].unique():
site_data = datasets[datasets['site'] == site]
maxdepth = site_data['z'].max()
for z in range(1, maxdepth + 1):
sample_data = site_data[site_data['z'] == z]
d_s_center = CalculateSauterDiameter(sample_data, 'd_center')
d_s_lower = CalculateSauterDiameter(sample_data, 'd_lower')
d_s_upper = CalculateSauterDiameter(sample_data, 'd_upper')
sauter_diameter = sauter_diameter.append({'site':site,
'z':z,
'd_s_center':d_s_center,
'd_s_lower':d_s_lower,
'd_s_upper':d_s_upper}, ignore_index=True)
sauter_diameter.to_csv('data/sauter_diameter.csv', index=False)
sauter_diameter.head(3)
d_s_center | d_s_lower | d_s_upper | site | z | |
---|---|---|---|---|---|
0 | 0.000011 | 0.000010 | 0.000011 | owens_lake_T8-W_P1 | 1.0 |
1 | 0.000018 | 0.000017 | 0.000019 | owens_lake_T8-W_P1 | 2.0 |
2 | 0.000005 | 0.000004 | 0.000005 | owens_lake_T8-W_P1 | 3.0 |