Learning Forum

[ehs3860]

Hi,

I want to simulate an ice accretion problem over a solid body.

The location at which ice appears and its height are defined via a DEFINE_GRID_MOTION macro;

The mesh displacement (according to the computed ice height) seems fine. However, since some of the regions over the body would address the presence of water over the ice (called Glaze Ice), I have no idea how to apply water presence on the ice and its temperature profile in these corresponding regions.

I tried to activate the VOF multiphase model and define the water as C_VOF(c,t) = 1.0 using a DEFINE_ADJUST macro. Still, I couldn't readily determine target cells in the interior domain perpendicular to the ice profile. Is my approach correct?

Any ideas would greatly be appreciated.

Kind Regards

[DrAmine]

What about just ignoring that water amount and consider it in your UDF too?

[ehs3860]

Thank you for your reply.
Do you mean this: "mesh displacement" = "ice-height" + "water-height"?
I'm afraid not; it is an aerodynamical problem over an airfoil, and ignoring the water film (or adding its height to that of the ice) could result in erroneous results.

[Rob]

How is the water film forming in the first place? Ie how do you get water onto the wing to then freeze?

[ehs3860]

Thank you, Rob, for your reply.
The icing model is computed using Myers's (2001) well-established model.
First, the supercooled water accumulation problem is solved using the Eulerian approach, and the corresponding differential equations are solved by UDS modules and the results are promising.
At the second stage, the previously obtained data are used in Myers's model to calculate the ice thickness as well as its temperature. Note that, depending on the physics of the flow, some regions may indicate only the ice (called Rime Ice), while other regions over the airfoil may address Glaze Ice formation, i.e., water over the ice. In either case, the airfoil profile is moved at the ice thickness magnitude computed at each computational cell inside the DIFINE_GRID_MOTION macro, with no errors. The problem I face is how to introduce the water presence to the solver.
Besides, in this regard, I've written a DEFINE_ADJUST macro for both the pressure and suction side of the airfoil which is as follows. The approach I used is:
loop over the airfoil face cells, to check whether the icing is the Glaze one (with water film)
If Glaze exists, loop over the computational cells (2D) and store the distance of any cells from the airfoil face and the angle between its vector from the airfoil's cell with the normal vector of the airfoil.
Check if the distance of the cell from the airfoil's face cell is within the computed water height over the airfoil, where the computed angle is minimum ( because not all cells are orthogonal to the airfoil wall; instead, find the best cell which is semi-perpendicular)
When the cell is within the water conditions, the temperature of the water film is calculated by Myers's and C_VOF(c,t) =1.0.
I am not sure if my approach is correct.
Another problem, the macro associated with the pressure side seems fine, but the one associated with the side caused the solver to crash with this error: "Node 1: Process 9924: Received signal SIGSEGV."
I don't know what the problem is.
SORRY FOR MY LONG REPLY...
DEFINE_ADJUST(water_film_ps,d)
{
real x[ND_ND];
real x_ps[ND_ND];
cell_t c, c0;
face_t f;
int ps_ID = 14; /*pressure side ID*/
Thread *ps = Lookup_Thread(d,ps_ID);
Thread *t0_ps = THREAD_T0(ps);
real NV_VEC(A_ps), A_ps_mag, NV_VEC(etha_ps), NV_VEC(coord_ps), dist_ps;
real ang_ps;
int clip_ID = 11; /*interior ID*/
Thread *clip = Lookup_Thread(d,clip_ID);
int clip_size = THREAD_N_ELEMENTS_INT(clip);
real (*arr_ps)[2] ;
arr_ps = (real (*)[2])malloc(2*clip_size*sizeof(real));

int loc_ps = 0;
int i = 0;

real T_f = 273.15; /*freezing temperature*/

/********/
begin_f_loop(f,ps)
{
F_CENTROID(x_ps,f,ps);
c0=F_C0(f,ps);
F_AREA(A_ps, f, ps);
A_ps_mag = NV_MAG(A_ps);
NV_VS(etha_ps,=,A_ps,/, -A_ps_mag);
/***check the GLAZE from the obtained data***/
if ((C_UDMI(c0,t0_ps,0)>0.0 && C_UDMI(c0,t0_ps,6)>0.0 && C_UDMI(c0,t0_ps,7)> 0.0) || (CURRENT_TIME-init_time) >= C_UDMI(c0,t0_ps,1))
{
begin_c_loop_int(c,clip)
{
C_CENTROID(x,c,clip);
NV_D(coord_ps, =, x[0]-x_ps[0], x[1]-x_ps[1], 0.0);
dist_ps =NV_MAG(coord_ps);
if (NV_DOT(coord_ps, etha_ps)>0.0)
{
ang_ps = fabs(acos(NV_DOT(coord_ps, etha_ps)/dist_ps));
}
else
{
ang_ps = 100.;
}

arr_ps[c][0] = ang_ps; /**assign computed anlges to the first column of array**/
arr_ps[c][1] = dist_ps; /**assign computed interior cell distance (from the airfoil face) to the second column of array**/
}
end_c_loop_int(c,clip)
/*** find minimum angle***/
real min_ang_ps = arr_ps[0][0];
for (i=0; i {
if (arr_ps[i][0] {
min_ang_ps = arr_ps[i][0];
loc_ps = i;
}
}
/**************/
if (arr_ps[loc_ps][1]<=C_UDMI(c0,t0_ps,6))
{
/***implemenation of the water presence*/
thread_loop_c(clip,d)
{
begin_c_loop_int(c,clip)
{
if (arr_ps[c][0] <= min_ang_ps && arr_ps[c][1]<=C_UDMI(c0,t0_ps,6) && C_UDMI(c0,t0_ps,6)>0.0)
{
C_VOF(c,clip) = 1.;
C_T(c,clip) = T_f + C_UDMI(c0,t0_ps,7) * arr_ps[c][1]; /** water film temp according to Myers's model**/
}
}
end_c_loop_int(c,clip)
}
}
}
}
end_f_loop(f,ps)

free(arr_ps);
}

[Rob]

Eulerian forming a film which then freezes makes sense. The trick may be to calculate the freezing bit and then move the mesh whilst simultaneously removing mass. That may not be overly easy, but should be possible. Personally, I'd be nudging for pointers on Fensap.....
Re the UDF, I have no idea, but it's traditional to include the headers which you're missing.