[NonStopEagle]

Hello everyone,
I have been trying to write a code for implementing maxwell first order slip boundary conitions at the diverging wall of a micro-nozzle. My code works well initially but eventually succumbs to sudden divergence. I believe i have figured out why this happens, the code is as follows :


#include "udf.h" /* must be at the beginning of every UDF */
/*#include "mem.h" /* must be at the beginning of this UDF */
//#include <math.h> /* must be at the beginning of this UDF */
/*#include "metric.h" /* must be at the beginning of this UDF */
#define k_B 1.3805e-23 /* Boltzman constant (j/K)*/
#define sigma 419e-12 /* molecular diameter*/
#define sigma_v 1 /* Tangential momentum accommodation coefficient*/
#define sigma_T 1 /*Thermal accommodation coefficient*/
#define pi 3.14159 /* pi number*/
#define sqrt_2 1.41421 /* sqrt(2) number */
#define len = 20e-5;
double relax = 0.1;
int counter = 0;
//#define relax 50e-3
DEFINE_PROFILE(wall_slip_velocity, t, index)
{
double T, P, density, MU, K, landa, a_m, du_da, dT_dx, u_f, DUY, DVY, DUX, DVX, u, v;
double A[ND_ND];
face_t f;
Thread* t0;
cell_t c0;
double DUR, DVR;


begin_f_loop(f, t)
{
t0 = THREAD_T0(t); /* adjacent cell thread to f */
c0 = F_C0(f, t); /* adjacent cell to f */
T = C_T(c0, t0); /* temperature of the cell */
P = C_P(c0, t0); /* peressure of the cell */
density = C_R(c0, t0); /* density of the cell*/
MU = C_MU_L(c0, t0); /* laminar viscosity of the cell */
u = C_U(c0, t0);
v = C_V(c0, t0);
/*DUR = C_U_RG(c0, t0)[1];
DVR = C_V_RG(c0, t0)[1];*/

K = C_K_L(c0, t0); /* thermal conductivity of the cell */
landa = (k_B * T) / (sqrt_2 * pi * sigma * sigma * P); /* mean free path line */

DUX = C_DUDX(c0, t0);
DVX = C_DVDX(c0, t0);
DUY = C_DUDY(c0, t0);
DVY = C_DVDY(c0, t0);
F_AREA(A, f, t);
a_m = NV_MAG(A); /* magnitude of the cell area vector */
du_da = NV_DOT(A, A) / a_m; /* n- component of the cell x-velocity reconstruction gradient(RG) vector */
dT_dx = 1; /* x-component of the cell temperature reconstruction gradient(RG) vector */
//F_PROFILE(f, t, index) = relax * (((2 - sigma_v) / sigma_v) * landa * (u / sqrt(u * u + v * v) * (((DUX + DVX) * (v / sqrt(u * u + v * v))) + ((DUY + DVY) * (u / sqrt(u * u + v * v)))) + 0.75 * (MU / (density * T)) * dT_dx);
//F_UDMI(f,t,1)= relax * (((2 - sigma_v) / sigma_v) * landa * (u / sqrt(u * u + v * v) * (((DUX + DVX) * (v / sqrt(u * u + v * v))) + ((DUY + DVY) * (u / sqrt(u * u + v * v)))) + 0.75 * (MU / (density * T)) * dT_dx);
F_UDMI(f, t, 0) = relax * (((2 - sigma_v) / sigma_v) * landa * (((DUX + DVX) * (v / sqrt(u * u + v * v))) - ((DUY + DVY) * (u / sqrt(u * u + v * v)))));
// F_PROFILE(f, t, index) = F_UDMI(f, t, 0);
//F_PROFILE(f, t, index) = F_UDSI(f, t, 0);
F_UDMI(f, t, 1) = landa;

}
end_f_loop(f, t)

}
DEFINE_PROFILE(temperature_jump, t, index)
{
face_t f;
cell_t c0;
Thread* t0;
double p, tem, g, lamda, num_den, temp, vis, temp1, x, i, mag, DT2;
double DT, prd, tc, tj, cp, a,u,v;
double QUANT;
g = 9.81;
a = 4.29866e-19;
begin_f_loop(f, t) // for each face: get face-id and thread
{
t0 = THREAD_T0(t); /* adjacent cell thread to f */
c0 = F_C0(f, t); /* adjacent cell to f */
p = C_P(c0, t0); // getting the pressure
tem = C_T(c0, t0); // getting the temperature
tc = C_K_L(c0, t0);
vis = C_MU_L(c0, t0);
cp = C_CP(c0, t0);
u = C_U(c0, t0);
v = C_V(c0, t0);
DT = C_T_G(c0, t0)[1];
DT2 = C_T_G(c0, t0)[0];
prd = cp * vis / (g * tc);

num_den = p / ((1.38064852e-23) * tem); // Calculation of the number density
temp1 = num_den * a;
lamda = 1 / (temp1);
QUANT = (((DT2) * (v / sqrt(u * u + v * v))) - ((DT) * (u / sqrt(u * u + v * v))));
tj = 380+relax *((((2 - sigma_v) / sigma_v) *1.14566*lamda * QUANT ));//0.15596 (380 + ((1.14566 * 1 * (lamda / prd) * ((0.9877 * DT)))) * 0.9877);
F_PROFILE(f, t, index) = tj;


}
end_f_loop(f, t)
}
DEFINE_EXECUTE_AT_END(define_relax)
{
if ((counter % 200) == 0)
{
relax = relax + 0.01;
printf("Relaxation : %f", relax);
}
counter = counter + 1;
}


In the code the gradients are calculated everytime the begin_f_loop() runs, and every change of the gradient values brings some error in the resulting slip velocity. This erred slip velocity produces further deviations in the gradients, until finally the solution becomes unstable and diverges.
The solution to this, as i believe is to run the simulation in no-slip conditions till convergence and then implement slip such that the gradient and other flow variables are calculated only once initially at all the faces of the slip wall and then to use those values to calculate the slip velocity.
And i need your help to accomplish this.
Any help regarding this will be highly appreciated.
thanks in advance.

[homay95]

Originally Posted by NonStopEagle

Hello everyone,
I have been trying to write a code for implementing maxwell first order slip boundary conitions at the diverging wall of a micro-nozzle. My code works well initially but eventually succumbs to sudden divergence. I believe i have figured out why this happens, the code is as follows :


#include "udf.h" /* must be at the beginning of every UDF */
/*#include "mem.h" /* must be at the beginning of this UDF */
//#include <math.h> /* must be at the beginning of this UDF */
/*#include "metric.h" /* must be at the beginning of this UDF */
#define k_B 1.3805e-23 /* Boltzman constant (j/K)*/
#define sigma 419e-12 /* molecular diameter*/
#define sigma_v 1 /* Tangential momentum accommodation coefficient*/
#define sigma_T 1 /*Thermal accommodation coefficient*/
#define pi 3.14159 /* pi number*/
#define sqrt_2 1.41421 /* sqrt(2) number */
#define len = 20e-5;
double relax = 0.1;
int counter = 0;
//#define relax 50e-3
DEFINE_PROFILE(wall_slip_velocity, t, index)
{
double T, P, density, MU, K, landa, a_m, du_da, dT_dx, u_f, DUY, DVY, DUX, DVX, u, v;
double A[ND_ND];
face_t f;
Thread* t0;
cell_t c0;
double DUR, DVR;


begin_f_loop(f, t)
{
t0 = THREAD_T0(t); /* adjacent cell thread to f */
c0 = F_C0(f, t); /* adjacent cell to f */
T = C_T(c0, t0); /* temperature of the cell */
P = C_P(c0, t0); /* peressure of the cell */
density = C_R(c0, t0); /* density of the cell*/
MU = C_MU_L(c0, t0); /* laminar viscosity of the cell */
u = C_U(c0, t0);
v = C_V(c0, t0);
/*DUR = C_U_RG(c0, t0)[1];
DVR = C_V_RG(c0, t0)[1];*/

K = C_K_L(c0, t0); /* thermal conductivity of the cell */
landa = (k_B * T) / (sqrt_2 * pi * sigma * sigma * P); /* mean free path line */

DUX = C_DUDX(c0, t0);
DVX = C_DVDX(c0, t0);
DUY = C_DUDY(c0, t0);
DVY = C_DVDY(c0, t0);
F_AREA(A, f, t);
a_m = NV_MAG(A); /* magnitude of the cell area vector */
du_da = NV_DOT(A, A) / a_m; /* n- component of the cell x-velocity reconstruction gradient(RG) vector */
dT_dx = 1; /* x-component of the cell temperature reconstruction gradient(RG) vector */
//F_PROFILE(f, t, index) = relax * (((2 - sigma_v) / sigma_v) * landa * (u / sqrt(u * u + v * v) * (((DUX + DVX) * (v / sqrt(u * u + v * v))) + ((DUY + DVY) * (u / sqrt(u * u + v * v)))) + 0.75 * (MU / (density * T)) * dT_dx);
//F_UDMI(f,t,1)= relax * (((2 - sigma_v) / sigma_v) * landa * (u / sqrt(u * u + v * v) * (((DUX + DVX) * (v / sqrt(u * u + v * v))) + ((DUY + DVY) * (u / sqrt(u * u + v * v)))) + 0.75 * (MU / (density * T)) * dT_dx);
F_UDMI(f, t, 0) = relax * (((2 - sigma_v) / sigma_v) * landa * (((DUX + DVX) * (v / sqrt(u * u + v * v))) - ((DUY + DVY) * (u / sqrt(u * u + v * v)))));
// F_PROFILE(f, t, index) = F_UDMI(f, t, 0);
//F_PROFILE(f, t, index) = F_UDSI(f, t, 0);
F_UDMI(f, t, 1) = landa;

}
end_f_loop(f, t)

}
DEFINE_PROFILE(temperature_jump, t, index)
{
face_t f;
cell_t c0;
Thread* t0;
double p, tem, g, lamda, num_den, temp, vis, temp1, x, i, mag, DT2;
double DT, prd, tc, tj, cp, a,u,v;
double QUANT;
g = 9.81;
a = 4.29866e-19;
begin_f_loop(f, t) // for each face: get face-id and thread
{
t0 = THREAD_T0(t); /* adjacent cell thread to f */
c0 = F_C0(f, t); /* adjacent cell to f */
p = C_P(c0, t0); // getting the pressure
tem = C_T(c0, t0); // getting the temperature
tc = C_K_L(c0, t0);
vis = C_MU_L(c0, t0);
cp = C_CP(c0, t0);
u = C_U(c0, t0);
v = C_V(c0, t0);
DT = C_T_G(c0, t0)[1];
DT2 = C_T_G(c0, t0)[0];
prd = cp * vis / (g * tc);

num_den = p / ((1.38064852e-23) * tem); // Calculation of the number density
temp1 = num_den * a;
lamda = 1 / (temp1);
QUANT = (((DT2) * (v / sqrt(u * u + v * v))) - ((DT) * (u / sqrt(u * u + v * v))));
tj = 380+relax *((((2 - sigma_v) / sigma_v) *1.14566*lamda * QUANT ));//0.15596 (380 + ((1.14566 * 1 * (lamda / prd) * ((0.9877 * DT)))) * 0.9877);
F_PROFILE(f, t, index) = tj;


}
end_f_loop(f, t)
}
DEFINE_EXECUTE_AT_END(define_relax)
{
if ((counter % 200) == 0)
{
relax = relax + 0.01;
printf("Relaxation : %f", relax);
}
counter = counter + 1;
}


In the code the gradients are calculated everytime the begin_f_loop() runs, and every change of the gradient values brings some error in the resulting slip velocity. This erred slip velocity produces further deviations in the gradients, until finally the solution becomes unstable and diverges.
The solution to this, as i believe is to run the simulation in no-slip conditions till convergence and then implement slip such that the gradient and other flow variables are calculated only once initially at all the faces of the slip wall and then to use those values to calculate the slip velocity.
And i need your help to accomplish this.
Any help regarding this will be highly appreciated.
thanks in advance.

hello! do you know how i could write udf cod for these functions below:
1. V= L(du/dy)
2. T)f = T)s + L(dT/dy)
T)f is fluid temperature and T)s is solid temperature and L is slip length
i would appreciate for any help. thanks

[vinerm]

You can use the code that you quoted in your last as well as current post. However, there is one common mistake in both of these codes; both use velocity to define slip boundary. That brings inconsistency and requires extra care from user end. The approach is to use flux profile, i.e., shear stress for velocity slip and heat flux for temperature slip. But if you are not well versed with UDFs, then, I'd recommend you use the code that you quoted in your last post. You can modify it as per your equations.

[homay95]

First of all thank you so much for your response, and yes I just started a project on michrochannels in ANSYS-FLUENT and I'm not professional in writing UDF codes. so you are saying that I can use this code as velocity slip and temperature jump for a solid-liquid interface?

[vinerm]

The code that you quoted in your first post. But do note that you will have to modify the equation being used since you want to use a different equation.

[homay95]

I understand but as I'm not so into programming, could you please be more specific? I mean which part of the code should be changed according to my equations?
Have you written a UDF for this (Maxwell's boundary condition) before?

[vinerm]

The UDF given at

UDF for slip boundary condition

is rather well documented. Final value applied is whatever is given on rhs of F_PROFILE. If you have C or C++ before, then modification would be easy. Even if you have not, then you can remove some of the lines within begin_f_loop{}end_f_loop but keep F_PROFILE line. Change the equation on its rhs with your equation. The gradient you need for your equation is already given in that code as strain rate magnitude.

[homay95]

okay thank you very much for your help. and there is just one more question: do you know how i could interpret more than one UDF to ANSYS-FLUENT?

[vinerm]

You can interpret only one file, however, no body stops you from putting any number of functions within the file. So, just put everything in single file.

[homay95]

okay. now I get it. and I'm really sorry for taking your time but as I said I'm starter in ANSYS-FLUENT, so I have so much questions. where should I apply this UDF for interface? I mean for the wall or for the shadow-wall between solid and liquid?

[vinerm]

Open the wall boundary and it will show you the cell zone it is associated with. You have to assign to the wall or its shadow, whichever is associated with the fluid zone.

[homay95]

So I should apply the UDF twice, once in the "Momentum-specified shear" as slip velocity and once in "Thermal-Temperature" as temperature jump, and after that should I leave the solid adjacent wall as "Coupled", is it right?

[vinerm]

Once you change it from coupled to temperature, it changes for both, i.e., for wall and its shadow.

[homay95]

Sorry but i didn't get it, after interpreting the UDF which includes both slip velocity and temperature jump in one file, where should I specifically assign this UDF file in the fluid zone. I mean for Momentum and Thermal seperately.

[vinerm]

Though you have single file, there are two functions, each with different name. You will see both names under all boundaries. Just select the name you want to apply.

[homay95]

okay then. thank you very very much for your help. I really owe you and sorry again for taking your time
best regards,
Homayoon

[homay95]

hello again. i've written a udf code for velocity slip boundary condition as follows:
#include "udf.h"
double dudy ,vel;
double current_vel;
DEFINE_PROFILE(slipvelocity, thread, position)
{
face_t f;
cell_t c;
Thread *tc;
double b = 0.000001;

begin_f_loop(f, thread)
{
c = F_C0(f,thread);
tc = THREAD_T0(thread);
dudy = C_U_G(c,tc)[1];
vel = b*dudy;
F_PROFILE(f, thread, position) = vel;
}
end_f_loop(f, thread)
}

is it correct? or does it need any modification?

[vinerm]

The code is syntactically correct, except that you have global variables. Define all variables within DEFINE_ macro and not outside it.

[homay95]

okay. thank you very much for your help. one more question is should i hook it as shear stress or moving wall velocity component in boundary conditions?

[vinerm]

Depends on what the UDF returns. Better approach is to write UDF for shear stress, however, your code returns velocity and should be used for velocity.