sandbox/alimare/phase_change_velocity.h
Phase change velocity calculation
The following set of functions modifies the phase change velocity field \mathbf{v}_{pc} = v_{pc} \mathbf{n}_{2\rightarrow 1} in interfacial cells and sets it to : \displaystyle \mathbf{v}_{pc} = \frac{1}{L_H}\left(\lambda_{1}\left.\nabla tr_{1}\right|_ {\Gamma} - \lambda_{2}\left.\nabla tr_{2}\right|_{\Gamma}\right) which is the Stefan condition, where L_H is the latent heat, tr_{1} and tr_{2} are scalars defined in phase 1 and phase 2 respectively and \lambda_i are the thermal diffusivities.
Anisotropy of the Gibbs-Thomson relation
We define 2 functions to take into account the anisotropy in the change velocity where we set the following condition : \displaystyle \epsilon = \overline{\epsilon} \left( 1. - \alpha \cos(n \theta + \theta_0)\right) where \theta is the angle between the interface and the x-axis, \alpha =0.5 and \theta_0 = 0 are hardcoded. This will be used in the Gibbs-Thomson formulation.
scalar Teq[];
#ifndef ANISO // if the user defines ANISO, then he wants to use standard
// anisotropic function, else he has to define its own fac1() function before
// the <<#include "phase_change_velocity.h">>
double fac1(coord n, double eps4){
if(eps4==0.)return 1.;
double sum = 0;
foreach_dimension()
sum += eps4*powf(n.x,4.); // taken from doi.org/10.1016/j.jcrysgro.2010.11.013
// simple fourfold anisotropy.
#if dimension ==2
double theta = atan2(n.y, n.x);
return 1.-15.*eps4*cos(4.*theta);
#else // DIMENSION == 3 && ANISO != 0
return (1.-3.*eps4+4*eps4*sum);
#endif
}
#endif // ifndef ANISO
For the gradient calculation on the interface, we need the temperature on the interface T_{\Gamma}, which is defined either by : \displaystyle T_{\Gamma} = T_{m} - \epsilon_{\kappa} * \kappa - \epsilon_{v} v_{pc} = T_m - \left(\overline{\epsilon_{\kappa}} * \kappa - \overline{\epsilon_{v}} v_{pc}\right)(1-0.5 cos(n \theta)) which is the classical Gibbs-Thomson equation.
double Temp_GT(Point point, double epsK, double epsV, vector vpc,
scalar curve, face vector fs, scalar cs, double eps4){
if(epsK==0. && epsV ==0.)return 0.;
double pref = 1.;
if(eps4 > 0.)pref = fac1(facet_normal( point, cs ,fs),eps4);
if(epsV==0.)return epsK*curve[]*pref;
else
{
double temp=0.;
foreach_dimension()
temp += sq(vpc.x[]);
return (epsK*curve[] - epsV*sqrt(temp))*pref;
}
}
Velocity in interfacial cells
This function modifies the vector v_pc[]
in interfacial cells using the Stefan condition. Note that this procedure is not sufficient to obtain a proper phase change velocity field used for a level-set, it must be complemented with a reconstruction procedure (see this page)
void phase_change_velocity_LS_embed (
scalar cs, face vector fs, scalar tr,
scalar tr2, double T_eq, vector v_pc,
double L_H, double lambda[2], double epsK,
double epsV, double eps4, scalar curve) {
We store in scalar Teq[]
the temperature on the interface.
#if Gibbs_Thomson
foreach(){
// here we suppose that Tm = 0
Teq[] = (interfacial(point, cs) ? T_eq + Temp_GT(point, epsK, epsV, v_pc,
curve, fs, cs, eps4) : nodata);
}
#else
foreach(){
Teq[] = (interfacial(point, cs) ? T_eq : nodata);
}
#endif
boundary({Teq});
restriction({Teq});
To calculate the gradient \left. \nabla tr \right|_{\Gamma}, we use the embed_gradient_face_x defined in embed that gives an accurate definition of the gradients with embedded boundaries.
vector gtr[], gtr2[];
boundary({tr});
foreach(){
foreach_dimension(){
gtr.x[] = 0.;
}
if(interfacial(point,cs)){
coord n, p;
embed_geometry(point, &p, &n);
double c = 0.;
double temp = Teq[];
double grad = dirichlet_gradient(point, tr, cs , n, p,
temp, &c);
For degenerate cases (stencil is too small), we add the term tr[]*c to the gradient.
foreach_dimension(){
gtr.x[] += grad*n.x+tr[]*c;
}
}
}
boundary((scalar*){gtr});
Now we need to change the volume fraction and face fractions to calculate the gradient in the other phase.
foreach(){
cs[] = 1.-cs[];
}
foreach_face()
fs.x[] = 1.-fs.x[];
boundary({cs});
restriction({cs});
boundary({tr2});
vector p_sauv[], n_sauv[];
foreach(){
foreach_dimension(){
gtr2.x[] = 0.;
}
if(interfacial(point,cs)){
coord n, p;
embed_geometry(point, &p, &n);
double c=0.;
double temp = Teq[];
double grad = dirichlet_gradient(point, tr2, cs , n, p,
temp, &c);
foreach_dimension(){
gtr2.x[] += grad*n.x+tr2[]*c;
}
}
}
boundary((scalar*){gtr2});
foreach(){
cs[] = 1.-cs[];
}
foreach_face()
fs.x[] = 1.-fs.x[];
boundary({cs,fs});
restriction({cs,fs});
With the the normal vector and the gradients of the tracers we can now compute the phase change velocity \mathbf{v}_{pc}.
foreach()
foreach_dimension()
v_pc.x[] = 0.;
foreach()
if(interfacial(point, cs))
foreach_dimension()
v_pc.x[] = (lambda[1]*gtr2.x[] - lambda[0]*gtr.x[])/L_H;
boundary((scalar *){v_pc});
}