Implementation of volume averaging in OpenFOAM


Nano-confined flows in presence of the wall form an inhomogeneous system, the potential part of the pressure tensor, `mathbf{P}^{U}(mathbf{r})` for such an inhomogeneous system based on Irving and Kirkwood formulation is

`mathbf{P}^U(mathbf{r}) = -frac{1}{2} langle sum_{i}^{N} sum_{ine j}^{N} frac{mathbf{r}_{ij}otimesmathbf{r}_{ij}}{r_{ij}} frac{d U(r_{ij})}{d r_{ij}} D(mathbf{r},mathbf{r}_{ij}) rangle`

where `mathbf{r}_{ij}=mathbf{r}_j – mathbf{r}_i` and `r_{ij} = |mathbf{r}_{ij}|` and `D(mathbf{r},mathbf{r}_{ij})` is the Dirac delta functional given as

`D(mathbf{r},mathbf{r}_{ij})= int_0^1ds delta(mathbf{r}_i – mathbf{r} +s mathbf{r}_{ij})`

The Dirac `delta(mathbf{r}_i – mathbf{r})` function gives the probability per unit volume of the i th molecule to be at `mathbf{r}` at time t. The variable `smathbf{r}_{ij}` `(0le sle1)` gives a location of the molecule from `mathbf{r}_i` to a location between the segments endpoints `mathbf{r}_i` and `mathbf{r}_j`. Thus the potential part manifests through molecular interaction between two particles and hence it is spatially non-local for domains of the order of the range of intermolecular forces. For an infinitesimal area the interacting pair whose line of centre meets at a point within the area contributes to the pressure tensor. Volume averaging method is one of the way to compute this spatially non-local potential part of the pressure tensor by localizing stress in an arbitrary shaped volume.

Implementation of volume averaging in OpenFOAM is quite easy due to its large base library, for example tensor and field operations are used to construct the potential part of the pressure tensor by an overloaded operator "*" which represents dyadic operator `otimes`. The molecules are represented by an object "atomisticMolecule", contribution to potential part of the pressure tensor from each pair (I,J) of molecule is given by mol_rf=fIJ*rIJ, where fIJ and rIJ are the force and separation vector for pair (IJ). To cut down the computation cost the pressure tensor calculation needs to be looped over the molecules, updating each bins with the contribution made by the corresponding interaction pair of molecules. Following is the code snippet of the volume averaging method for slab bins as shown in the figure

  
//looping for all molecules I
for (int i = 0; i < Total_mols; i++)
{
molI = Occupancy[i];
// population of molecules in each bin
get_bin_pop(bin_pop);
for (int j = 0; j < Total_mols; j++)
{
//looping for all molecules J > I
if(j > i)
{
molJ = Occupancy[j];
evaluate_force(molI, molJ, fIJ);
vector rIJ = molI->position() - molJ->position();
// potential part as dyadic product of fIJ and rIJ
mol_rf = fIJ*rIJ;
scalar mol_delz = molJ->position().z()-molI->position().z();
get_bin_molIJ(k_top,k_bot,bm_topz,bm_botz);
for(int k=k_top+1; k < k_bot; k++)
{
bin_rf[k]+= mol_rf*bin_delz/mol_delz;
}
//end points of the segment rIJ in k_top and k_bot
bin_rf[k_top]+=mol_rf*bm_topz/mol_delz;
bin_rf[k_bot]+=mol_rf*bm_botz/mol_delz;
}
}
}
//bin averaging
for(int k=k_top; k < k_bot+1; k++)
{
bin_rf[k] = bin_rf[k]/bin_pop[k];
}

The potential part of the pressure tensor contribution to the kth bin is updated in tensorfield bin_rf and bin_topz, bin_botz are the segments in the partially intercepted first and the last bin as shown in the figure.