About the basic idea about PNM solving process.

[xu-kai-xu]

According to what I have learned in articles and in OpenPNM, I found PNM solves this equation for pressure distribution:

$$ Q = \sum\limits_{i = 1}^N {{q_i}}  = \sum\limits_{i = 1}^N {\sum\limits_{j = 1}^n {{g_{ij}}({P_j} - {P_i})} } = 0 $$

Is this correct?

And I found in OpenPNM, all the data would be transferred in to a matrix Ax=b. I can understand A denotes the hydraulic conductances in each pore location. And x denotes pressure at each pore location. But I don't know what b means. Is this understanding right?

Thanks for your help.

For example, this article introduced the PNM calculation details:

Gostick, J.T., Ioannidis, M.A., Fowler, M.W., et al. Pore network modeling of fibrous gas diffusion layers for polymer electrolyte membrane > fuel cells. Journal of Power Sources, 2007, 173(1): 277-290.

[jgostick]

b is the boundary conditions. In the equation above b is all 0's, but this would lead to a trivial solution. On your boundary pores you write that b = R as the rate of material entering or exiting a pore. For internal pores this balances to 0, but for boundary pores this is non-zero. If you have a 'value' boundary condition, then the row in the A matrix is chaged by be x_i=1, so b=1.

[xu-kai-xu]

Thanks, @jgostick, then I make a demo for this, in which I set a 2D PNM with 3*4=12 pores, with pore diameter being 50 micron, and throat diameter being 25 microns, lattice spacing 100 microns. Below is the pore and throat connected.

The hydraulic conductance is 1.80689199e-13 because all throats are the same.

I set the inlet pores are [3, 7, 11], with pressure 20 Pa, outlet pores are [0, 4, 8], with pressure 10 Pa. Then I got the matrix A and b.

flow.set_value_BC(pores=inlets, values=20)
flow.set_value_BC(pores=outlets, values=10)

A = 
  (1, 2)	-1.8068919869130543e-13
  (5, 6)	-1.8068919869130543e-13
  (9, 10)	-1.8068919869130543e-13
  (1, 5)	-1.8068919869130536e-13
  (2, 6)	-1.8068919869130536e-13
  (5, 9)	-1.8068919869130543e-13
  (6, 10)	-1.8068919869130543e-13
  (2, 1)	-1.8068919869130543e-13
  (6, 5)	-1.8068919869130543e-13
  (10, 9)	-1.8068919869130543e-13
  (5, 1)	-1.8068919869130536e-13
  (6, 2)	-1.8068919869130536e-13
  (9, 5)	-1.8068919869130543e-13
  (10, 6)	-1.8068919869130543e-13
  (0, 0)	5.119527296253653e-13
  (1, 1)	5.420675960739162e-13
  (2, 2)	5.420675960739162e-13
  (3, 3)	5.119527296253653e-13
  (4, 4)	5.119527296253653e-13
  (5, 5)	7.227567947652215e-13
  (6, 6)	7.227567947652216e-13
  (7, 7)	5.119527296253653e-13
  (8, 8)	5.119527296253653e-13
  (9, 9)	5.420675960739162e-13
  (10, 10)	5.420675960739163e-13
  (11, 11)	5.119527296253653e-13

b = 
   5.11952730e-12, 1.80689199e-12, 3.61378397e-12, 1.02390546e-11,
   5.11952730e-12, 1.80689199e-12, 3.61378397e-12, 1.02390546e-11,
   5.11952730e-12, 1.80689199e-12, 3.61378397e-12, 1.02390546e-11

I can not find the relation between b and the boundary conditions. May you explain that in more details?

Flow equation in OpenPNM.ipynb.txt

I also upload the ipynb file. And wish this could be useful as a demo for new beginners.

[jgostick]

Hi @xu-kai-xu , I'm sorry, my previous answer was a rushed and incomplete...I didn't realize you'd be digging in so deeply! When you set a value BC, we do as I said above, but also two other tricks for numerical efficiency.

Firstly, we scale the BC values so that the A matrix remains well scaled (i.e. we don't want 1.0 and 1e-15 in the same matrix). To do this, we compute the average value of the diagonal, and then multiply the A and b values by this average, so the numbers you see in b are not exactly the value BCs you entered, but are multiplied by this factor. So instead of x_i=BC_i, the equation we enter into the system is f*x_i=f*BC_i where f is mean(A.diagonal()).

Secondly, all pores that are internal but connected to a boundary pore get an additional adjustment. We subtract the BC value from A and add it to b. This may not be necessary, and I'm going to open an issue to discuss whether we should remove it.

BTW, your efforts to understand the inner workings are very impressive. I always tell my students to do a whole calculation by hand to see how it all works, which you are doing. If you want to make your job easier, I suggest simplifying even further. Here is a very basic script that will help:

pn = op.network.Cubic([3, 3, 1])
a = op.phases.GenericPhase(network=pn)
a['throat.hydraulic_conductance'] = 1.0
sf = op.algorithms.StokesFlow(network=pn, phase=a)
sf.set_value_BC(pores=[0, 1, 2], values=1.0)
sf.set_value_BC(pores=[6, 7, 8], values=0.0)
print(sf.b)  # will be all zeros since we haven't actually set the BC yet
print(sf.A.todense())  # will be easily readable since the g_h values are just 1.0's
diag = A.diagonal()  # Let's get a copy of the diagonal, which is used to compute f
f = np.mean(diag)
print(diag)  # This has not been altered yet
sf._apply_BCs()  # A private method that get's called when you run sf
print(sf.b)  # Will now have some values in it
diag = A.diagonal()
print(diag)  # This has been altered by the application of the value BCs
print(sf.b / f)  # This should give the b values as I described previously, showing the BCs