Pore/throat widening due to dissolution and diffusion?

[nltthinh]

Hi, it's me, the newbie, again. I tried to search in discussion but didn't find anything.

I'm trying to simulate the dissolution and diffusion of minerals in a piece of rock (leaching process). In reality, an aggressive solution (i.e. water, acid) dissolves the solid mineral from the pore wall, then dissolved ions diffuse out of the pores. This leaching process enlarges the pores/throats that leads to higher porosity, and eventually higher permeability.

To simulate such a thing I need to:

• Provide a value of concentration in each pore/throat. This concentration represents the ions that are formed by the dissolution of the solid mineral.
• Those ions will diffuse out of the pores, and further dissolution occurs to replenish this outflow.
• The dissolution will widen pores/throats.

Please give me an idea on how to proceed. Many thanks in advance.

[jgostick]

What you are trying to do is possible with openpnm but is a bit advanced. Basically, you'll want to do an advection-diffusion simulation with reactions in the pores. This will get you started. Assuming that your reaction is relatively slow, then you can do 'quasi' transient, meaning you solve the steady-state problem, find the rate of dissolution (from the reaction term in each pore), multiply this rate by some time step (1 day, 1 week or whatever), then compute how much mass was lost. THEN you'll have to get creative, and find how much bigger each pore got based on how much mass it lost, then recompute their sizes accordingly. Then you regenerate all the size dependent models and repeat the above process. We don't currently have a built in couple of algorithms with changing geometry, so this is why you need to do it by hand.

[nltthinh]

Hi it's me again,

Since I ignore the advection part, I decide to use Fickian to simulate the diffusion of ions out of the pores. The "leaching" process (aka dissolution + diffusion) is usually assumed to be diffusion controlled since dissolution is fast. What I'm trying to do is:

• define a function to calculate concentration of ions in the pores based on a parameter called CaSi, which is the mole ratio of Calcium over Silicon of the solid. This function reproduces a curve of [Ca] in pores vs Ca/Si ratio of solid that I took from literature.
• deploy fickian diffusion. After some time step, CaSi will be recalculated from the change in [Ca], and new value of [Ca] is renewed from the mentioned function. pore.concentration will be saved in an array for later access.
• the pore enlargement part I've not found a way yet.

Below is my code, please replace the pore.diameter with either webull or lognormal.
For now I'm stuck at the tfd.run() step, as it gives the error "Found NaNs in A matrix, possibly caused by NaNs in throat.diffusive_size_factors, throat.hydraulic_size_factors. The issue might get resolved if you call regenerate_models on the following object(s): geo_01, phase_01"

ws = op.Workspace()
ws.settings['loglevel'] = 40
np.set_printoptions(precision=4)

# create pore network
r = 5e-2                    #in m
h = 2.5e-2                  #in m
cyl = op.network.DelaunayVoronoiDual(num_points=(1000), shape=[r,h])

geo = op.geometry.GenericGeometry(network=cyl, pores=cyl.Ps, throats=cyl.Ts)
geo.add_model(propname='pore.seed',
              model=op.models.geometry.pore_size.random,
              seed=0)
geo['pore.diameter'] = D
geo.add_model(propname='pore.cross_sectional_area',
                       model=op.models.geometry.pore_cross_sectional_area.sphere,
                       pore_diameter='pore.diameter')
geo.add_model(propname='pore.volume',
                       model=op.models.geometry.pore_volume.sphere,
                       pore_diameter='pore.diameter')
geo.add_model(propname='throat.max_size',
                       model=op.models.misc.from_neighbor_pores,
                       mode='min',
                       prop='pore.diameter')
geo.add_model(propname='throat.diameter',
                       model=op.models.misc.scaled,
                       factor=0.5,
                       prop='throat.max_size')
geo.add_model(propname='throat.length',
                        model=op.models.geometry.throat_length.spheres_and_cylinders,
                        pore_diameter='pore.diameter',
                        throat_diameter='throat.diameter')
geo.add_model(propname='throat.cross_sectional_area',
                        model=op.models.geometry.throat_cross_sectional_area.cylinder,
                        throat_diameter='throat.diameter')
geo.add_model(propname='throat.volume',
                        model=op.models.geometry.throat_volume.cylinder,
                        throat_diameter='throat.diameter',
                        throat_length='throat.length')
geo.add_model(propname='throat.diffusive_size_factors',
                        model=op.models.geometry.diffusive_size_factors.spheres_and_cylinders,
                        pore_diameter="pore.diameter",
                        throat_diameter="throat.diameter")
geo.add_model(propname='throat.hydraulic_size_factors',
                        model=op.models.geometry.hydraulic_size_factors.spheres_and_cylinders,
                        pore_diameter="pore.diameter",
                        throat_diameter="throat.diameter")

# create phases
water = op.phases.GenericPhase(network=cyl)
water['pore.temperature'] = 293
water.add_model(propname='pore.viscosity', model=op.models.phases.viscosity.water)

# check permeability
P1 = 6.05e5
P2 = 0
water.add_model(propname='throat.hydraulic_conductance',
                model=op.models.physics.hydraulic_conductance.classic_hagen_poiseuille)
sf = op.algorithms.StokesFlow(network=cyl)
sf.setup(phase=water)
sf.set_value_BC(pores=cyl.pores('top'), values=P1)
sf.set_value_BC(pores=cyl.pores('bottom'), values=P2)
sf.run()

Q = sf.rate(pores=cyl.pores('top'))     #m3/s
A = np.pi*((r)**2)                      #m2
L = h                                   #m
mu = 0.001                              #Pa.s
dP = P1-P2                              #Pa
K = Q*mu*L/A/dP
print(f'Permeability is {K}m2 at dP = {dP}Pa or {dP/100000}bar and Q = {Q}m3/s or {Q*6e10}uL/min')

water.update(sf.results()) 

# define liquid concentratrion vs Ca/Si
def c_Ca_liquid(CaSi):
    global c_Ca_liquid_max
    if CaSi > 1.895:
        c_Ca_liquid_max = 20.7e-3
    elif (CaSi <= 1.895) and (CaSi > 1.2):
        c_Ca_liquid_max = 23.62*CaSi-24.07
    elif (CaSi <= 1.2) and (CaSi > 0.67):
        c_Ca_liquid_max = 1.00757+0.00953*np.exp(4.85486*CaSi)
    elif (CaSi <= 0.67) and (CaSi > 0.123):
        c_Ca_liquid_max = 1.27
    elif (CaSi <= 0.123) and (CaSi >= 0):
        c_Ca_liquid_max = 1.31707-1.31446*np.exp(-26.2786*CaSi)
    return c_Ca_liquid_max

# initialize Ca/Si
density_ini = 2.3e6   #g/m3
CaO = 0.6192          #%CaO in CEM I
SiO2 = 0.1947         #%SiO2 in CEM I
V_pores = sum(geo['throat.volume']) + sum(geo['pore.volume']) 
V_sample = A*L
wc = 0.5
nCasolid = density_ini*V_sample/(1+wc)*CaO/56.04        # mol
nSisolid = density_ini*V_sample/(1+wc)*SiO2/60.08       # mol

# diffusivity, pore concentration
water['pore.diffusivity'] = water['throat.diffusivity'] = 8e-10     # m2/s
water.add_model(propname='throat.diffusive_conductance',
                     model=op.models.physics.diffusive_conductance.generic_diffusive,
                     pore_diffusivity="pore.diffusivity", 
                     throat_diffusivity="throat.diffusivity", 
                     size_factors="throat.diffusive_size_factors")
water.regenerate_models()
tfd = op.algorithms.TransientFickianDiffusion(network=cyl, phase=water)
tfd.setup(t_scheme='cranknicolson', t_final=1, t_output=1, t_step=1, t_tolerance=1e-12)
tfd.set_value_BC(pores=cyl.pores('boundary'), values=0)
c_Ca_pore = []
c_Ca_liquid_max = 0
for i in range(1,2): 
    CaSi = nCasolid/nSisolid
    c_Ca_liquid_max = c_Ca_liquid(CaSi)
    tfd.set_IC(values=c_Ca_liquid_max)     
    tfd.run()
    c_Ca_pore.append(tfd['pore.concentration'])
    nCasolid = nCasolid-(tfd['pore.concentration@0']-tfd['pore.concentration'])*V_pores

[jgostick]

The code above looks ok, so my best guess if that you've got some odd numbers in your geometry that is causing the conductance models to break. For instance if your pores are large and they overlap by too much, then the conductance can become nan. You need to track down the original source of the nans. I'm guessing if you print(geo) you'll see that the size factor arrays have some nans. To fix this you need to be more careful about assign pore sizes. There is a model in openpnm.models.geometry.pore_size which computes the largerst sphere you can put in each pore without overlapping a neigbhor, so many use this as a check?