Learning Forum

[ilya.nearakcheev]

I am trying to model the Marangoni effect on a free surface where heating occurs. In my problem using two materials, but implicitly the third is air. So I want to use the Marangoni condition on the surface:

Hovewer in Ansys 2023R1 Fluent does not support non-"single-valued" expressions to account for the distribution of different phases on the surface, that is, the coefficient must be geometrically the same on the edge. The most interesting thing is that the assignment of conditions on the shear stress vector is supported by expressions, but for some reason I get an error about "unsuccessful interpolation in the nodes of the face", therefore there is simply zero by default, as can be seen in the solution in the form of a wall shear stress equal to 0.
So I took a stab at creating a UDF that works, but the problem is that I also use internal sources via DEFINE_SOURCE that regulates the speed in cells that have not yet melted .And here is the strangest thing, that only the Marangoni stress condition built into Fluent gives correct results, and the velocities are also "suppressed" by the sources. In the case DEFINE_PROFILE is basically applie boundary conditions on the velocity, without taking into account the sources included in the momentum equation (my guess), unlike the built-in method. In this regard, I provide my code and ask for help with how exactly you implement your conditions so that, as you write, "This shear stress is then applied to the momentum equation" (apparently directly).

#include "vars_opt_conv_rad.h"

#include "udf.h"

#include "sg.h"

#include "sg_mphase.h"

#include "mem.h"

#include "metric.h"

#include "math.h"

 

#define dTdx 0

#define dTdz 1

#define VOF1 2

#define LVOF1 3

#define LVOF2 4

 

DEFINE_EXECUTE_ON_LOADING(on_loading, libudf)

{

Set_User_Memory_Name(dTdx, "dTdx");

Set_User_Memory_Name(dTdz, "dTdz");

Set_User_Memory_Name(VOF1, "VOF1");

Set_User_Memory_Name(LVOF1, "LVOF1");

Set_User_Memory_Name(LVOF2, "LVOF2");

}

void compute_temperature_derivatives(Domain *domain)

{

    MD_Alloc_Storage_Vars(domain, SV_T_RG, SV_T_G, SV_NULL);

 

    Scalar_Reconstruction(domain, SV_T, -1, SV_T_RG, NULL);

 

    Scalar_Derivatives(domain, SV_T, -1, SV_T_G, SV_T_RG, NULL);

}

 

DEFINE_ADJUST(compute_grad_and_phase, domain)

{

    Thread *ct, *phase_thread;

    cell_t c;

real *Tgr;

real T;

 

    compute_temperature_derivatives(domain);

 

    thread_loop_c(ct, domain)

    {

        begin_c_loop(c, ct)

        {

 

            Tgr = C_T_RG(c, ct); 

            C_UDMI(c, ct, dTdx) = Tgr[0]; 

            C_UDMI(c, ct, dTdz) = Tgr[2];  

 

            phase_thread = THREAD_SUB_THREAD(ct, 0);  

            if (phase_thread)

                C_UDMI(c, ct, VOF1) = C_VOF(c, phase_thread);

 

real Tm = C_UDMI(c, ct, VOF1) * Tm_phase1 + (1.0 - C_UDMI(c, ct, VOF1)) * Tm_phase2;

 

T = C_T(c, ct);

 

if (T <= Ts1) {

C_UDMI(c, ct, LVOF1) = 0.0; 

} else if (T >= Tl1) {

C_UDMI(c, ct, LVOF1) = 1.0; 

} else {

C_UDMI(c, ct, LVOF1) = (T - Ts1) / (Tl1 - Ts1); 

}

 

if (T <= Ts2) {

C_UDMI(c, ct, LVOF2) = 0.0;

} else if (T >= Tl2) {

C_UDMI(c, ct, LVOF2) = 1.0;

} else {

C_UDMI(c, ct, LVOF2) = (T - Ts2) / (Tl2 - Ts2); 

}

        }

        end_c_loop(c, ct)

    }

}

 

DEFINE_PROFILE(marangoni_shear_x, t, i)

{

    face_t f;

    cell_t c;

    Thread *ct;

 

    begin_f_loop(f, t)

    {

        c = F_C0(f, t);

        ct = F_C0_THREAD(f, t);

 

        real vof1 = C_UDMI(c, ct, VOF1);  

real f1 = C_UDMI(c, ct, LVOF1);

real f2 = C_UDMI(c, ct, LVOF2);

        real gradTx = C_UDMI(c, ct, dTdx);

        real T = C_T(c, ct);               

 

        real Tm = vof1 * Tm_phase1 + (1.0 - vof1) * Tm_phase2;

 

        real dSigma_dT = D_SIGMA_DT_1 * vof1 + D_SIGMA_DT_2 * (1.0 - vof1);

real f_liquid = vof1 * f1 + (1.0 - vof1) * f2;

 

        real sigmoid = f_liquid / (1.0 + exp((Tm - T) / 10));

 

        F_PROFILE(f, t, i) = - dSigma_dT * gradTx * sigmoid;

    }

    end_f_loop(f, t)

}

 

DEFINE_PROFILE(marangoni_shear_z, t, i)

{

    face_t f;

    cell_t c;

    Thread *ct;

    

    begin_f_loop(f, t)

    {

        c = F_C0(f, t);

        ct = F_C0_THREAD(f, t);

 

        real vof1 = C_UDMI(c, ct, VOF1);   

real f1 = C_UDMI(c, ct, LVOF1);

real f2 = C_UDMI(c, ct, LVOF2);

        real gradTz = C_UDMI(c, ct, dTdz);

        real T = C_T(c, ct);

 

        real Tm = vof1 * Tm_phase1 + (1.0 - vof1) * Tm_phase2;

 

        real dSigma_dT = D_SIGMA_DT_1 * vof1 + D_SIGMA_DT_2 * (1.0 - vof1);

real f_liquid = vof1 * f1 + (1.0 - vof1) * f2;

 

        real sigmoid = f_liquid / (1.0 + exp((Tm - T) / 10));

 

        F_PROFILE(f, t, i) = - dSigma_dT * gradTz * sigmoid;

    }

    end_f_loop(f, t)

}

 

 

/* Brake force X */

DEFINE_SOURCE(s_bfx, c, t, dS, eqn)

{

    real src, T, vx0;

    real velocity_multiplier, brake_factor;

    real f1, f2;

    Thread *primary_thread;

 

    primary_thread = THREAD_SUB_THREAD(t, 0);

    real phase1_frac = C_VOF(c, primary_thread);

    T = C_T(c, t);

    real Tm = phase1_frac * Tm_phase1 + (1.0 - phase1_frac) * Tm_phase2;

    

    f1 = C_UDMI(c, t, LVOF1);

    f2 = C_UDMI(c, t, LVOF2);

    

    brake_factor = brakeForceFactor * (

        phase1_frac * (1.0 - f1) +  

        (1.0 - phase1_frac) * (1.0 - f2)

    );

    

    vx0 = 0;  

    velocity_multiplier = -C_R(c, t) * brake_factor / (1.0 + exp((T - Tm) / 3.0));

    

    src = velocity_multiplier * (C_U(c, t) - vx0);

    dS[eqn] = velocity_multiplier;

    

    return src;

}

 

/* Brake force Y */

DEFINE_SOURCE(s_bfy, c, t, dS, eqn)

{

    real src, T, vy0;

    real velocity_multiplier, brake_factor;

real f1, f2;

    Thread *primary_thread;

 

primary_thread = THREAD_SUB_THREAD(t, 0);

real phase1_frac = C_VOF(c, primary_thread);

T = C_T(c, t);

real Tm = phase1_frac * Tm_phase1 + (1.0 - phase1_frac) * Tm_phase2;

 

f1 = C_UDMI(c, t, LVOF1);

f2 = C_UDMI(c, t, LVOF2);

    brake_factor = brakeForceFactor * (

        phase1_frac * (1.0 - f1) + 

        (1.0 - phase1_frac) * (1.0 - f2) 

    );

 

    vy0 = 0;  

 

    velocity_multiplier = -1.0 * C_R(c, t) * brake_factor / (1.0 + exp((T - Tm) / 3.0));

 

    src = velocity_multiplier * (C_V(c, t) - vy0);

    dS[eqn] = velocity_multiplier;

 

    return src;

}

 

/* Brake force Z */

DEFINE_SOURCE(s_bfz, c, t, dS, eqn)

{

    real src, T, vz0;

    real velocity_multiplier, brake_factor;

    real f1, f2;

    Thread *primary_thread;

 

    primary_thread = THREAD_SUB_THREAD(t, 0);

    real phase1_frac = C_VOF(c, primary_thread);

    T = C_T(c, t);

    real Tm = phase1_frac * Tm_phase1 + (1.0 - phase1_frac) * Tm_phase2;

    

    f1 = C_UDMI(c, t, LVOF1);

    f2 = C_UDMI(c, t, LVOF2);

    

    brake_factor = brakeForceFactor * (

        phase1_frac * (1.0 - f1) +  

        (1.0 - phase1_frac) * (1.0 - f2)

    );

    

    vz0 = 0;  

    velocity_multiplier = -C_R(c, t) * brake_factor / (1.0 + exp((T - Tm) / 3.0));

    

    src = velocity_multiplier * (C_W(c, t) - vz0);

    dS[eqn] = velocity_multiplier;

    

    return src;

}

[Rob]

Staff won't be able to comment on the coding. 

For the source terms, you may find they're calculated sequentially so freezing the flow is done and then shear is added: you may need to do some experimentation. Where are you applying the DEFINE_PROFILE? Surface tension is DEFINE_PROPERTY for some models, and not sure on the phase interaction panel. 

[ilya.nearakcheev]

I am trying to model the interaction of a separate phase - two metals with essentially a quasi-phase of air. That is, we are not talking about the interaction between the volume fractions of two metals, which are in the problem, where I've already set DEFINE_PROPERTY in interphase interaction. Therefore, instead of introducing a third phase (air), I set a boundary condition on the free surface, on which, among other things, heating occurs. And within the framework of this physical model, I need to write DEFINE_PROFILE for "Specific Shear Boundary Conditions" where it's possible set UDF. 

[Rob]

OK, so the wall is there instead of the air. And you're then fixing the cell velocity adjacent to the wall? 

[ilya.nearakcheev]

 

Actually yes, the speed is through shear stress vector. It is on the upper face that I am trying to get the Marangoni effect, as originally intended in your built-in boundary condition “Marangoni Stress”

 

[Rob]

Which may conflict with the wall function if you're then trying to freeze the neighbouring cell. I assume wall adhesion is off? 

[ilya.nearakcheev]

If you mean adhesion in VOF, then yes, and there’s no other kind either.

In general, when I look at the results and check the velocity contours on that surface, the velocities are quite high and incorrect — unlike the built-in "Marangoni stress" boundary condition. However, if I disable Boundary Values on the velocity contour, the result matches the built-in one exactly (I tested this using the same d sigma/dT coefficient).

I even tried to further suppress the velocities near the neighboring cells adjacent to the face, but it didn’t help. That’s why I thought that the boundary condition is applied independently and just influences the solution, but in essence, the surface velocities are “physically” incorrect.

[Rob]

So, the surface values in your UDF agree with the built in model, and it's the cell zone effects that you're adding that are causing the issue? 

[ilya.nearakcheev]

 

No, I’m using suppressing sources via DEFINE_SOURCE both with the built-in “Marangoni Stress” boundary condition and with my own DEFINE_PROFILE, and the surface velocities are correct only with the built-in one.

But as I said, its use is limited to a single-value expression. That’s why I started all this in the first place — I don’t understand where the error is, or perhaps Fluent itself cannot properly include a DEFINE_PROFILE into the momentum equation, if I understood the description of the built-in “Marangoni Stress” correctly from the earlier discussion.