sandbox/fintzin/Rising-Suspenion/Drop.h
Declaration of the structures gathering all the properties linked to a droplet.
- The property of the drop are indexed and just “a” in the code.
typedef struct {physical properties
double vol; // volume
double S; // surface
double p; // pressure
double W2; // internal fluctuation
coord Y; // position of the center of mass
coord U; // averaged velocity
coord Uc; // centered velocity
coord Fh; // drag force
coord D_nbr; // distance to the nearest neigboring particles
coord Per; // Indicator to identify teh peeriodicity of a drop
tens P; // moment of momentum tensor defined as $\int \bm{ru} dV$
tens G; // shape tensor defined as $\int \bm{rr} dV$
double Gnorm; // norm of the inertia tensor
double Pnorm; // norm of the moment of momentum tensor
tens WW; // tensor product of the inner fluctuation tensor around Uc- @brief Think about how to compute the first moments of the surface
- forces body forces and hydrodynamic forces
- You might add the int of the stress here
- so that you can compute the sum of each moments.
- What about the centered grad of velocity ?
tens Mh; // First hydrodinamical moment.
tens Sig; // Sigma stress tesnor
tens Ms; // Moment of the surface tension
double ed; // energy dissipation- bilan on NRJ
Other properties linked to the identification
double ct; // contact time with the nearest neigbor
int tagmin; // tag of the nearest droplet
int tag;
int tagmax;
int realtag;
int si;
int adj[100]; // maximum index of scalar feilds
} Drop;assign a value to detect if a drop cross a bounday
void isperio(scalar s, scalar tag ,int nb,int tagi, Drop n[] ){
coord per1[nb],per2[nb];
for (int i = 0; i < nb; i++)
foreach_dimension()
per1[i].x = 0,per2[i].x = 0;
foreach(reduction(+:per1[:nb]) reduction(+:per2[:nb]))
if (s[] != nodata && dv() > 0. && s[]>=EPS) {
int i=tag[]-1;
coord pos = POS;
foreach_dimension(){
if(pos.x>(L0/2-0.2*D))
per1[i].x = 1;
else if(pos.x<-(L0/2-0.2*D))
per2[i].x = 1;
}
}
for (int i = 0; i < nb; i++)
foreach_dimension()
n[tagi+i].Per.x = (per1[i].x && per2[i].x);
}compute the center of mass and average velocity of the drop
void pos_vol_and_vel(scalar s, scalar tag ,int nb,int tagi,Drop n[]){
double volume[nb];
coord posG[nb],v[nb];
for (int i = 0; i < nb; i++){
volume[i]=0;
posG[i]=v[i]=Coord_nul;
}
foreach(reduction(+:posG[:nb]) reduction(+:v[:nb]) reduction(+:volume[:nb]))
if (s[] != nodata && dv() > 0. && s[]>=EPS){
int i=tag[]-1;
volume[i] += dv()*s[];
coord pos = POS;
pos = POS_perio(pos,n[tagi+i].Per);
foreach_dimension(){
posG[i].x += dv()*s[]*pos.x;
v[i].x += dv()*s[]*u.x[];
}
}
for(int i=0;i<nb;i++){
n[tagi+i].vol = volume[i];
n[tagi+i].Y = volume[i] ? div_coord(posG[i],volume[i]) : Coord_nul;
n[tagi+i].U = volume[i] ? div_coord(v[i],volume[i]) : Coord_nul;
}
}store the value of the velocity feild at the center of mass location for each drops
void points_vel(scalar s, scalar tag ,int nb,int tagi,Drop n[]){
for(int i=0;i<nb;i++){
coord velc_loc={0},velc_global = {0};
#if dimension == 2
double a = n[i].Y.x ,b=n[i].Y.y;
foreach_dimension()
velc_loc.x = interpolate(u.x,a,b) < 100*normL2_coord(n[i].U) ? interpolate(u.x,a,b) : 0. ;
#elif dimension == 3
double a = n[i].Y.x ,b=n[i].Y.y, c = n[i].Y.z;
foreach_dimension()
velc_loc.x = interpolate(u.x,a,b,c) < 20*normL2_coord(n[i].U) ? interpolate(u.x,a,b,c) : 0. ;
#endif
#if _MPI
MPI_Allreduce(&velc_loc.x, &velc_global.x, dimension, MPI_DOUBLE, MPI_SUM,MPI_COMM_WORLD);
#else
velc_global = velc_loc;
#endif
n[tagi+i].Uc = velc_global;
}
}compute the inner velocity w_i dV turbulent tensor \int w_i w_j dv and \int w^c_i dV
void compute_WW_W2(scalar s, scalar tag ,int nb,int tagi,Drop n[]){
tens WW[nb];
double W2[nb];
for (int i = 0; i < nb; i++) {
WW[i] = tens_nul;
W2[i] = 0;
}
foreach (reduction (+:WW[:nb]) reduction(+:W2[:nb]))
if (s[] != nodata && dv() > 0. && s[]>=EPS) {
int i=tag[]-1;
einstein_sum(i,j){
WW[i].i.j += (u.i[] - n[tagi+i].U.i)*(u.j[] - n[tagi+i].U.j)*dv()*s[];
W2[i] += (u.i[] - n[tagi+i].U.i)*(u.j[] - n[tagi+i].U.j)*dv()*s[];
}
}
for(int i=0;i<nb;i++){
n[tagi+i].WW = WW[i];
n[tagi+i].W2 = W2[i];
}
}compute the moment of momentum of the drop P = r_i u_j dV
void compute_G_and_P(scalar s, scalar tag ,int nb,int tagi,Drop n[]){
// need to check the periodicity
tens G[nb],P[nb];
for (int i = 0; i < nb; i++)
G[i]=P[i]=tens_nul;
foreach (reduction (+:G[:nb]) reduction (+:P[:nb]))
if (s[] != nodata && dv() > 0. && s[]>=EPS) {
int i=tag[]-1;
coord Pos = POS;
coord r = dist_perio_coord(Pos,n[tagi+i].Y);
einstein_sum(i,j){
G[i].i.j += r.i * r.j * dv() * s[];
P[i].i.j += r.i * u.j[] * dv() * s[];
}
}
for(int i=0;i<nb;i++){
n[tagi+i].G = G[i];
n[tagi+i].P = P[i];
n[tagi+i].Gnorm = tens_norm(G[i]);
n[tagi+i].Pnorm = tens_norm(P[i]);
}
}Compute the integral of the newtonian stress tensor
void compute_Mh(scalar s, scalar tag ,int nb,int tagi,Drop n[]){
tens Mh[nb],Sig[nb];
double ed[nb],pressure[nb];
for (int i = 0; i < nb; i++){
ed[i] = 0;
pressure[i] = 0;
Mh[i] = tens_nul;
Sig[i] = tens_nul;
}
foreach (reduction (+:Sig[:nb]) reduction (+:Mh[:nb])
reduction (+:ed[:nb]) reduction(+:pressure[:nb]))
if (s[] != nodata && dv() > 0. && s[]>=EPS) {
int i = tag[]-1;
pressure[i] += p[]*s[]*dv();
tens gradU,S,rdivS,rdivp;
coord gradp,LapU;
coord Pos = POS;
coord r = dist_perio_coord(Pos,n[tagi+i].Y);
vec_grad(u,gradU);
vec_lap(u,LapU);
scalar_grad(p,gradp);
// need to think about what tensor is needed
einstein_sum(i,j){
S.i.j = (gradU.i.j+gradU.j.i)*mu1;
rdivp.i.j = r.i * gradp.j;
rdivS.i.j = r.i * LapU.j * mu1;
Sig[i].i.j += S.i.j * s[] * dv();
Mh[i].i.j += (S.i.j - rdivp.i.j + rdivS.i.j)*s[]*dv();
ed[i] += S.i.j*gradU.i.j*s[]*dv();
}
}
for(int i=0; i<nb; i++){
n[tagi+i].Mh = Mh[i];
n[tagi+i].Sig = Sig[i];
n[tagi+i].ed = ed[i];
n[tagi+i].p = pressure[i];
}
}compute the integral of the surface tension force and its moment
void compute_S_and_Ms(scalar s, scalar tag ,int nb,int tagi,Drop n[]){
double S[nb];
tens Ms[nb];
for (int i = 0; i < nb; i++){
Ms[i] = tens_nul;
S[i] = 0;
}
foreach(reduction(+:S[:nb]) reduction(+:Ms[:nb])){
int i = tag[]-1;
if(s[] != nodata && dv() > 0. && s[] > EPS && s[] < 1. - EPS){
int i = tag[]-1;
coord normal = mycs(point, s), p;
double alpha = plane_alpha(s[], normal);
double dS = pow(Delta, dimension - 1) * plane_area_center(normal, alpha, &p);
S[i] += dS;
}
if(i < 0){
foreach_neighbor(1)
i = max(i,tag[]-1);
}
if(i>=0){
coord Pos = POS;
coord r = dist_perio_coord(Pos,n[tagi+i].Y);
coord acc = {a.x[1],
a.y[0,1],
a.z[0,0,1]};
coord grad = {s[1] != s[-1],
s[0,1] != s[0,-1],
s[0,0,1] != s[0,0,-1]};
einstein_sum(i,j){
Ms[i].i.j += (grad.j && fm.j[] > 0) ?
r.i * (acc.j+a.j[])/2. * (f[]*(rho1-rho2)+rho2) * dv() :0.;
}
}
}
for(int i=0; i<nb; i++){
n[tagi+i].S = S[i];
n[tagi+i].Ms = Ms[i];
}
}
void Compute_drops (scalar s, scalar tag ,int nb,int tagi,Drop n[]){
isperio(s,tag,nb,tagi,n);
pos_vol_and_vel(s,tag,nb,tagi,n);
for(int i=0;i<nb;i++)
foreach_dimension()
if (n[tagi+i].Y.x>L0/2.)
n[tagi+i].Y.x -= L0; // set the pos inside the domain
points_vel(s,tag,nb,tagi,n);second order term calculation appearing in the * moment of momentum equation
compute_G_and_P(s,tag,nb,tagi,n);
compute_WW_W2(s,tag,nb,tagi,n);
compute_Mh(s,tag,nb,tagi,n);
compute_S_and_Ms(s,tag,nb,tagi,n);
}funciton that link every drops from one step to the next by the variable tag
void assign_tags_from_last_step (Drop datat0[],Drop datat1[],int tagi){
int tags_avail[tagi];
for(int k = 0;k<tagi;k++){
datat1[k].tag = -1;
tags_avail[k]=0;
}
for (int j = 0; j < datat0[0].tagmax; j++)
{
coord np = add_coord(datat0[j].Y,mult_coord(datat0[j].U,dtprint)); // euleur schem
#if dimension == 2
double r = sqrt(datat0[j].vol/M_PI)*0.5; //
#else
double r = pow(datat0[j].vol/M_PI,1./3.)*3./4.; //
#endif
bool probleme_of_assignement = true;
for (int k = j; k < (tagi+j); k++)
{
int idx = k % tagi;
coord np1 = datat1[idx].Y;
double dist = dist_perio(np,np1);
if(dist<r*0.5){
datat1[idx].tag = datat0[j].tag;
probleme_of_assignement = false;
break;
}
}
if(probleme_of_assignement){
fprintf(stdout,"Pbug : X %g, Y %g, Z %g\n",np.x,np.y,np.z);
fprintf(stdout,"Error of assignement at t = %g for bubble tag = %d tm %d j %d vol %g \n",
t,datat0[j].tag,datat0[0].tagmax,j,datat0[j].vol);
}
}to reassgin tags when coalescence and breakup event occure
int ta = 0;
for (int k = 0; k < tagi; k++){
bool avail = true;
for (int j = 0; j < tagi; j++)
if (k == datat1[j].tag)
avail = false;
if(avail) {
tags_avail[ta] = k;
ta++;
}
if(datat1[k].tag>=tagi) datat1[k].tag=-1;
}
int new_tag=0;
for (int k = 0; k < tagi; k++){
if(datat1[k].tag == -1 && new_tag<=ta){
if(tags_avail[new_tag] >= tagi) new_tag++;
datat1[k].tag = tags_avail[new_tag];
new_tag++;
}
}
}print the distance between each droplets and store the tag of the closest * dorplet of a given droplet. besides it store the contact time
void print_pair_dist(Drop data0[],Drop data[],int tagi,int step, float t){
Drop data_orderred[tagi];
for (int j = 0; j < tagi; j++)
{
data[j].ct = 0;
foreach_dimension()
data[j].D_nbr.x = Ls*100;
data_orderred[data[j].tag]=data[j];
}
foreach_dimension(){
Dist_x = fopen(name_x,"a");
fprintf(Dist_x,"%d,%g",step,t);
}
for (int k = 0; k < tagi; k++)
{
coord npi = data_orderred[k].Y;
for (int j = k+1; j<tagi;j++){
coord npj = data_orderred[j].Y;
coord distcoord = dist_perio_coord(npj,npi);
double dist = normL2_coord(distcoord);
foreach_dimension()
fprintf(Dist_x,",%g",distcoord.x);
if(dist < normL2_coord(data_orderred[k].D_nbr)){
data_orderred[k].D_nbr = distcoord;
data_orderred[k].tagmin = j;
}
if(dist < normL2_coord(data_orderred[j].D_nbr)){
data_orderred[j].D_nbr = mult_coord(distcoord,-1);
data_orderred[j].tagmin = k;
}
}
for (int j = 0; j < tagi; j++)
{
if(data[j].tag == k) {
data[j].D_nbr = data_orderred[k].D_nbr;
data[j].tagmin = data_orderred[k].tagmin;
}
}
}
foreach_dimension(){
fprintf(Dist_x,"\n");
fclose(Dist_x);
}
#if COAL
Drop data_orderred0[data0[0].tagmax];
// compute the contact time
for (int j = 0; j < data0[0].tagmax; j++)
data_orderred0[data0[j].tag]=data0[j];
for(int k = 0; k<tagi;k++){
int tag = data[k].tag;
int tm = data[k].tagmin;
double Dia = sqrt(data_orderred[tm].vol/ M_PI )+sqrt(data[k].vol/ M_PI );
if((normL2_coord(data[k].D_nbr) <= Dia) && (data[k].adj[data_orderred[tm].si] == true)){
if(tag >= data0[0].tagmax) data[k].ct = dtprint;
if(tag < data0[0].tagmax) data[k].ct = data_orderred0[tag].ct + dtprint;
}
}
#endif
}Function that print the property of the droplets in fDrop.csv and * print the pair dist in Dist_x.csv
void print_Drop(Drop data0[],Drop data1[],int tagi,int step,double time){
fDrop = fopen("fDrop.csv","a");
for(int j=0;j<tagi;j++){
fprintf(fDrop,"%d,%g,%d",step,time,data1[0].tagmax);
fprintf(fDrop,",%g",data1[j].vol);
fprintf(fDrop,",%g",data1[j].S);
fprintf(fDrop,",%g",data1[j].ed);
fprintf(fDrop,",%g",data1[j].Gnorm);
fprintf(fDrop,",%g",data1[j].Pnorm);
fprintf(fDrop,",%g",data1[j].p);
fprintf(fDrop,",%g",data1[j].W2);
einstein_sum(i,j){
fprintf(fDrop,",%g",data1[j].Y.i);
fprintf(fDrop,",%g",data1[j].U.i);
fprintf(fDrop,",%g",data1[j].Uc.i);
fprintf(fDrop,",%g",data1[j].D_nbr.i);
fprintf(fDrop,",%g",data1[j].Fh.i);
fprintf(fDrop,",%g",data1[j].G.i.j);
fprintf(fDrop,",%g",data1[j].P.i.j);
fprintf(fDrop,",%g",data1[j].WW.i.j);
fprintf(fDrop,",%g",data1[j].Sig.i.j);
fprintf(fDrop,",%g",data1[j].Mh.i.j);
fprintf(fDrop,",%g",data1[j].Ms.i.j);
}
fprintf(fDrop,",%g",data1[j].ct);
fprintf(fDrop,",%d",data1[j].tagmin);
fprintf(fDrop,",%d",data1[j].tag);
fprintf(fDrop,"\n");
}
fclose(fDrop);
}
#if COAL
void neighboring_scalar(scalar s, scalar tag ,int nb,int tagi,Drop n[]){
int taggg[nb];
int nvar = datasize/sizeof(double), len = nb*nvar, adj[len],tc[len];
for (int i = 0; i < len; i++) adj[i] = false;
for (int i = 0; i < nb; i++) taggg[i]=tc[i]=0;
foreach(reduction(max:taggg[:nb]) reduction(max:tc[:nb])
reduction(max:adj[:nvar]))
if (s[] != nodata && dv() > 0. && s[]>=EPS) {
int i=tag[]-1;
taggg[i] = tag[];
foreach_neighbor()
for(scalar c in interfaces)
if(c.i != s.i && c[] > EPS)
adj[i*c.i + c.i] = true;
}
for (int i = 0; i < nb; i++){
n[tagi+i].realtag = taggg[i];
for(scalar c in interfaces){
n[tagi+i].adj[c.i] = adj[i*c.i + c.i];
}
}
}coalescence function. * It is activated only if COAL = 1
void coalescence(Drop data[],int tagi){
Drop DAT[tagi];
for (int j = 0; j < tagi; j++) DAT[data[j].tag]=data[j];
scalar b[];
for(scalar c in interfaces){
foreach(){
b[] = c[]>EPS;}
tag (b);
int si[tagi],realtag[tagi];
for (int i = 0; i < tagi; i++) si[i]=realtag[i]=-1;
for(int j = 0 ;j < tagi ; j++){
if(DAT[j].ct != 0.){
fprintf(stdout,"time %g the tag %d %d ct %g\n",t,j,DAT[j].tag,DAT[j].ct);
Drop b1 = DAT[j];
Drop b2 = DAT[DAT[j].tagmin];
double d_eq = 2*sqrt(b2.vol/M_PI)*sqrt(b1.vol/M_PI) / (sqrt(b2.vol/M_PI) + sqrt(b1.vol/M_PI));
double V_0 = normL2_coord(diff_coord(b1.U,b2.U));
double td = rho_f*sq(d_eq)*V_0/(8*sig);
double We = rho_f*d_eq*sq(V_0)/(2*sig);
fprintf(stdout,"coalescence time td = %g We = %g deq %g vel %g\n",td,We,d_eq,V_0);
if (DAT[j].ct > td && DAT[j].si != DAT[DAT[j].tagmin].si && c.i == DAT[j].si) {
// for(int l = 0 ;l < tagi ; l++){
// int tm = DAT[j].tag;
// if(DAT[l].tag != DAT[j].tag && DAT[l].tag != DAT[tm].tag){
// if(DAT[l].adj[DAT[j].si] == 1 && DAT[l].si == DAT[tm].si)
// if(DAT[l].tagmin == DAT[j].tag){
// fprintf(stdout,"l too close too [");
// for(scalar s in interfaces) fprintf(stdout,"(%d,%d),",DAT[l].adj[s.i],s.i);
// fprintf(stdout,"]\n");
// }
// }else{
realtag[j] = DAT[j].realtag;
si[j] = DAT[DAT[j].tagmin].si;
DAT[j].si = DAT[DAT[j].tagmin].si;
fprintf(stdout,"add to caaol realt %d si objc %d\n",si[j],realtag[j]);
// }
// }
}
}
}
#if _MPI
MPI_Bcast(si,tagi,MPI_INT,0,MPI_COMM_WORLD);
MPI_Bcast(realtag,tagi,MPI_INT,0,MPI_COMM_WORLD);
#endif
for (int i = 0; i < tagi; i++)
if(si[i] != -1)
foreach()
if(b[] == realtag[i]){
scalar t = {si[i]};
t[] = c[];
c[] = 0.;
}
scalar * list = list_copy ({c});
for (int i = 0; i < tagi; i++)
if(si[i] != -1)
list = list_add(list,(scalar){si[i]});
boundary (list);
free (list);
}
}
#endif