sandbox/tavares/contact_embed.h
Contact angles on an embedded body
This file is used to impose contact angles on arbitrary embedded bodies described using [EMBED] (embed.h).
//#include "fractions.h"
#include "embed.h"
#include "geometrical_tools.h"
#define Pi 3.14159265358979323846
extern scalar tag, triple_cell, f1, Alpha;
We compute the normal corresponding to the angle we want to impose on the contact cell knowing the initial volume fraction in that cell.
coord normal_contact (coord ns, coord nf1, double angle){
double det = ns.x*nf1.y - ns.y*nf1.x; // 2D cross product
coord n;
if ( det > 0){
n.x = ns.x*cos(angle)-ns.y*sin(angle); // new fluid normal with the right contact angle
n.y = ns.x*sin(angle)+ns.y*cos(angle);
}
else{
n.x = ns.x*cos(angle)+ns.y*sin(angle);
n.y = -ns.x*sin(angle)+ns.y*cos(angle);
}
return n;
}
void triple_cell_dectection(double angle, scalar f, scalar cs, vector ns, vector nf1, vector nf){
scalar alpha1[],alpha2[];
foreach(){
tag[] = triple_cell[] = 0;
int n=0;
foreach_neighbor(){
if(cs[]<1. && cs[]>0. && f[]<1. && f[]>0. ){
We locate the exact cell containing the triple point
coord b, ntemp, ntemp1, ntemp2, p_vof[2], p_embed[2];
vof geometry
ntemp1 = interface_normal (point, f); // normal to the fluid interface
foreach_dimension() nf1.x[] = ntemp1.x ;
//alpha1[] = line_alpha(f[], ntemp1);
alpha1[] = line_alpha(f[],(coord){nf1.x[], nf1.y[]});
embed solid geometry
embed_geometry (point, &b, &ntemp2);//ntemp = facet_normal(point, cs, fs);
//embed_geometry (point, &b, (coord*){ns.x[], ns.y[]});
//embed_geometry (point, &b, &ns[]);
foreach_dimension() ns.x[] = -ntemp2.x ;
alpha2[] = line_alpha(cs[], ntemp2);
Does VOF interface intersects the embed solid in this cell ?
if ( facets (ntemp1, alpha1[], p_vof)==2 && facets (ntemp2, alpha2[], p_embed)==2 ){
if (Intersect(p_vof[0],p_vof[1],p_embed[0],p_embed[1])) {
Then the triple point lies in here
triple_cell[]=1; n++;
double piangle = (angle*Pi)/180. ;
//foreach_dimension() ntemp2.x = ns.x[] ;
ntemp = normal_contact((coord){ns.x[], ns.y[]}, (coord){nf1.x[], nf1.y[]}, piangle);
foreach_dimension() nf.x[] = ntemp.x;
Alpha[] = line_alpha(f[], ntemp); // fictive alpha is stored
}
else nf.x[] = nf.y[] = ns.x[] = ns.y[] = nodata;
}
}
}
if(n && cs[]==0.) tag[] =1.;
}
}
void fraction_reconstruction(scalar f, vector nf){
foreach(){
if(tag[]){
int stencil=2;
for (int i1 = -stencil; i1 <= stencil; i1++){
for (int j1 = -stencil; j1 <= stencil; j1++){
if(triple_cell[i1,j1] ){
coord a, b;
a.x = -0.5 - i1;
a.y = -0.5 - j1;
b.x = 0.5 - i1;
b.y = 0.5 - j1;
f1[] = rectangle_fraction ((coord){nf.x[i1,j1], nf.y[i1,j1]}, Alpha[i1,j1], a, b);
}
}
}
}
}
boundary({f1});
foreach(){
if (tag[]){
f[]=f1[];
}
}
boundary({f});
}