Calculation of interface area interFoam / VoF
 

Peace Foamers,


I'm currently implementing an evaporation model and for that, I need the interface Area between my two phases alpha1 and alpha2.


Does anyone know how to do that? I assumed that one of the interface reconstruction methods would deliver a function that could be called or the class surfaceInterpolation but I havent found anything suitable there. (btw. I'm not greatly familiar with PLIC, isoAdvector, etc.)


My first thought of an apporach looked like that:
1.) Identifying interface cells (e.g. alpha = 0.5)
2.) Then calculate the volume occupied by alpha for each of these interface cells (mesh.V()*alpha1)

3.) calculate normal vector with nHat = grad(alpha)/mag(grad(alpha)) to get information about the orientation of that approximated interface plane
______________ Now the part where I dont know how to proceed
Here I thought about an algorithm that is starting from a suitable point. I try to explain it in a simple case: 2D Square. Starting from a corner then follow the diagonal and calculate a orthogonal line where both ends of this line touch the sides of that square. This will be continued until the areas (equivalent to the volume in 2.)) will match. The resulting length(area) is what can be used for the calculation.


Has someone an idea how to realize that - or even better: Is there something which gives me the interface area between both phases without a complicate implementation?


Kind regards and a thank you in advance!

Found a solution?
 

Hey überschwupper,


basically I need a pretty similar procedure to detect collisions between two flow fronts of alpha1... (track interface (get element attributes to current interface elements etc,). You already found a good way?


There must be some function in interFoam to get the interface elements (eg. based on the gradient of alpha?).



Tank you in advance!
Julian

Hey, there is a function that already tracks interface elements (coming from the interfaceProperties class inhereted by incompressibleTwoPhaseMixture - nearInterface()) but anyway.. This interface area calculation is still open for me.

Hi uberschwupper,
I have a similar doubt. I wnat to find the distance between cell center and interface. Do you have any idea about that?
Thanks in advance

mag(fvc::grad(alpha1))*mesh().V()

Hi Santiago
mag(fvc::grad(alpha1))*mesh().V()it has the dimensions of area. fvc::grad(alpha1) has dimensions 1/length and mesh.V() has dimensions Length*Length*Length. Do you know how to find the distace from cell center to interface.
Thanks in advance

Denote the interface by \Gamma, and follow the usual FOAM conventions. If P is the owner cell, and N is of the neighbour then the coordinate to the free-surface from the centroid of an interface cell is equal to:

X_\Gamma = X_P + \zita times d_f

where

\zita = \frac{\alpha_P -0.5}{\alpha_P - \alpha_N}

Hi santiago,
I am new to foam can I have a detailed explaination? I don’t know what alpha_n alpha_p and d_f are.
Thanks and regards

X means coordinate; and d_f = X_N - X_P, \alpha is the indicator function (alpha1, alpha.water, or whatnot).

Hi santiago,
I mean P and N indicate alpha1 at owner and neighbour cells?

Can you please post equation of zeta in mathematical form instead of code?
Thanks and regards

Yes.

\zita = (\alpha_P - 0.5)/(\alpha_P - \alpha_N)

P and N indices are for current cell and neighbour cell right?

for the 2nd time, yes.

Thank you santiago I will try this out .

Hi santiago ,
Can you please tell me what does zita physically means? Can I get the physical significance of zita?
I implemented code in this way but I am getting large value of center to interface distance . value is around 500. Can you look into it and suggest any changes? Or if you have any code sample please give me.


const unallocLabelList& neighbour = mesh.neighbour();
const unallocLabelList& owner = mesh.owner();
const volVectorField& C = mesh.C();
const volScalarField& T = alpha1().db().lookupObject<volScalarField>("T");
//const surfaceScalarField& deltas = mesh.deltaCoeffs();
//const dimensionedScalar len("len",dimLength,1/mag(deltas));
volScalarField zita ("zita", 0*limalpha1);
vector delta(0,0,0);
const dimensionedScalar len("len",dimLength,0);
const dimensionedScalar smallValue("smallValue",dimLength,1e-5);
const dimensionedScalar secondValue("secondValue",dimensionSet(1,-2,-1,0,0,0,0),0.0 );
//const volScalarField thirdValue("thirdValue",1e-6 * limalpha1);
forAll(owner, facei)
{
zita = (limalpha1[owner[facei]] - 0.5)/(max((limalpha1[owner[facei]] - limalpha1[neighbour[facei]]),value));
vector delta = C[neighbour[facei]] - C[owner[facei]];
const dimensionedScalar len("len",dimLength,mag(delta));
//Info << "zita: " << max(zita).value() << " ND" <<endl;
//Info << "len: " << len.value() << " m" <<endl;
}

return Pair<tmp<volScalarField>>
(
(limalpha1*lambda1_+ (1-limalpha1)*lambda2_) * ((T-Tsat_)/(max((mag(C)+(zita*len)),smallValue)*h_lv_)),limal pha1 * 0.0 * secondValue
);
}

Thanks in advance

There is none. The formula comes from writing a linear interpolation operator to the interface and zita is the weight of said interpolation, between an owner and a neighbour cell that are crossed by the interface. This formula is not valid anymore in non-interface cells, hence the error you are getting.

If you need the distance to the interface from ALL CELLS then you not only need to determine the distance to the interface centroids, but you need to solve a Poisson equation for coordinate quantities wrt the interface location. The faceMeshWave, and meshWave class may help you with that. You may get some inspiration from some wall & turbulence models where the distance from the wall is needed.

Hi santiago,
I don’t need distance from every cell . Just from the cell center to the cell containing interface. Can you attach any code sample if possible? or please tell me how to get X_p and alpha values of owner and neighbor cells in openfoam.
Thanks in advance

But you are doing it in your code! The only problem is that you are not correctly identifying your interface cells. An interface cell is one where the product (\alpha_P - 0.5)*(\alpha_N - 0.5) < 0.

Hi santiago,
So I should use for loop instead of forAll .
I am not getting idea how to implement it with for loop. Because If I implement it with for loop I will have 2 return statements. One for non interface cells and other is for interface cells.
Can you attach code snip if you have it?
Thanks in advance