[BW1211]

Hello everyone. Recently，I am using the DPM model to simulate the ash deposition in the boiler. I used the critical viscosity model in udf to describe ash deposition process, and the deposition layer was divided into three layers: particle inner layer, sintered layer and molten layer. It is worth noting that when the surface temperature reaches a certain value, the particle layer can transform to a sintered layer and the sintered layer can transform to a molten layer. I want to get the deposition thickness of each layer through simulation.(This is the background)
My problem is that when I use the udf below, the deposition thickness of the sintered layer is negative (UDM3), which is unreasonable.I hope you can help me check the udf (I have checked it many times, but I haven't found the problem). Thank you for your help.
DEFINE_DPM_BC(capture_bounce_bc, p, t, f, f_normal, dim)
{

real p_loc[ND_ND],p_temperature, p_diameter;
real plocy,pvel;
real A[ND_ND],area,flow,flux;
real viscosity,visc_exponent;
real p_rep_mass,deposit_mass;

Thread *t0;
cell_t c0;
c0=F_C0(f,t);
t0=F_C0_THREAD(f,t);

//Calculate the mass represented by the impacting particle

p_diameter=P_DIAM(p); //Get Particle Diameter from Fluent
p_rep_mass=RR_MASS_FLOW*TIMESTEP_SIZE*(RR_DIAM_MAX-RR_DIAM_MIN)*(RR_SF/p_diameter)*pow((p_diameter/RR_DBAR),RR_SF)*pow(EXP,-1*(pow((p_diameter/RR_DBAR),RR_SF)))/RR_NUM_PARTICLES;


//Determine viscosity of ash particle

p_temperature=P_T(p); //Get particle temperature from fluent
visc_exponent=14788/(p_temperature-T_Shift)-10.931; //Calculate Ash viscosity
viscosity=(p_temperature-T_Shift)*pow(10,visc_exponent); // Finish Viscosity calculation



//Assume deposited mass is inversely proportional to viscosity when particle viscosity is higher than critical viscosity

if(viscosity>critical_viscosity)
{
deposit_mass=p_rep_mass*critical_viscosity/viscosity;
}

else
{
deposit_mass=p_rep_mass;
}



//Calculate mass deposited per unit face area

F_AREA(A,f,t);
area=NV_MAG(A); //Get face area from fluent
flux=deposit_mass/area;

F_UDMI(f,t,0)+=deposit_mass;
/*C_UDMI(c0,t0,0)=F_UDMI(f,t,0);*/

//Add thickness to slag layer if the face temperature is higher than T_Slag, then kill the particle and return

if(F_T(f,t)>T_Slag)
{

F_UDMI(f,t,4)+=flux/RHO_Slag;
/*C_UDMI(c0,t0,4)=F_UDMI(f,t,4);*/
return(PATH_ABORT);

}

//Add thickness to sintered layer if the face temperature is higher than T_Sinter, then kill particle and return

if(F_T(f,t)>T_Sinter)
{
F_UDMI(f,t,3)+=flux/RHO_Sinter;
/*C_UDMI(c0,t0,3)=F_UDMI(f,t,3);*/
return(PATH_ABORT);
}

//Add thickness to particulate layer if the surface temperature is lower than T_Sinter, then kill particle and return

F_UDMI(f,t,2)+=flux/RHO_Particulate;
/*C_UDMI(c0,t0,2)=F_UDMI(f,t,2);*/
return(PATH_ABORT);

}



DEFINE_DPM_INJECTION_INIT(Injection_Init,I)
{

int i;
real current_time=0;
Thread *t;
Domain *d;
face_t f;
real R_Part,R_Sint,R_Slag, R_Total, L_Part, L_Sint, L_Slag, T_SintSlag_ifc,
T_PartSint_ifc, T_Surface, Delta_L_Sint,
Delta_L_Melt,Y_Sint,Y_Part;
d=Get_Domain(1);




//Check Sintering and Melting

thread_loop_f(t,d)
{

if(NNULLP(THREAD_STORAGE(t,SV_UDM_I)))
{

begin_f_loop(f,t);

{

L_Part=F_UDMI(f,t,2);
L_Sint=F_UDMI(f,t,3);
L_Slag=F_UDMI(f,t,4);



R_Part=L_Part/K_Particulate;
R_Sint=L_Sint/K_Sinter;
R_Slag=L_Slag/K_Slag;


R_Total=R_Part+R_Sint+R_Slag+R_BASE;
T_Surface=F_T(f,t);


cell_t c0;
Thread *t0;
c0=F_C0(f,t);
t0=F_C0_THREAD(f,t);



//Calculate interface temperatures between ash layers

T_PartSint_ifc=T_Cool+(T_Surface-T_Cool)*(R_Part+R_BASE)/R_Total;
T_SintSlag_ifc=T_Cool+(T_Surface-T_Cool)*(R_BASE+R_Part+R_Sint)/R_Total;

//Ckeck Melting and melt appropriate amount of sintered layer

if(T_SintSlag_ifc>T_Slag)
{

Y_Sint=R_Total*K_Sinter*(T_Slag-T_Cool)/(T_Surface-T_Cool)-K_Sinter*(R_Part+R_BASE);
Delta_L_Melt=L_Sint-Y_Sint;
F_UDMI(f,t,3)=Y_Sint;
F_UDMI(f,t,4)=L_Slag+Delta_L_Melt*RHO_Slag/RHO_Sinter;

C_UDMI(c0,t0,3)=F_UDMI(f,t,3);
C_UDMI(c0,t0,4)=F_UDMI(f,t,4);
}

//Check sintering and sinter appropriate amount of particulate layer

if(T_PartSint_ifc>T_Sinter);
{

Y_Part=R_Total*K_Particulate*(T_Sinter-T_Cool)/(T_Surface-T_Cool)-R_BASE*K_Particulate;
Delta_L_Sint=L_Part-Y_Part;
F_UDMI(f,t,2)=Y_Part;
F_UDMI(f,t,3)=L_Sint+Delta_L_Sint*RHO_Sinter/RHO_Particulate;

C_UDMI(c0,t0,2)=F_UDMI(f,t,2);
C_UDMI(c0,t0,3)=F_UDMI(f,t,3);


}


}

end_f_loop(f,t)

}

}




//Increment simulation time, print to screen and store current simulation time at UDM NO.5

thread_loop_f(t,d)
{

if(NNULLP(THREAD_STORAGE(t,SV_UDM_I)))

{

begin_f_loop(f,t)

{
cell_t c0;
Thread *t0;
c0=F_C0(f,t);
t0=F_C0_THREAD(f,t);

F_UDMI(f,t,5)+=TIMESTEP_SIZE;
/*C_UDMI(c0,t0,5)=F_UDMI(f,t,5);*/


}

end_f_loop(f,t)

}
}

}

[pakk]

Your UDF is big. Can you remove the parts that are not relevant for this problem?

[BW1211]

Thank you for your reply. The udf can compile well and the main problem is that I use UDMI(2,3,4) as the calculation variable for the subsequent slag layer conversion, and the value of UDMI3 is negative. I hope you can help to see if there is a problem between these F_UDMIs. Thanks a lot. The simplified udf looks like this：
DEFINE_DPM_BC(capture_bounce_bc, p, t, f, f_normal, dim)
{
deposit_mass=p_rep_mass;
if(viscosity>critical_viscosity)
{
deposit_mass=p_rep_mass*critical_viscosity/viscosity;
}
//Calculate mass deposited per unit face area
F_AREA(A,f,t);
area=NV_MAG(A); //Get face area from fluent
flux=deposit_mass/area;

F_UDMI(f,t,0)+=deposit_mass;
C_UDMI(c0,t0,0)=F_UDMI(f,t,0);

//Add thickness to slag layer if the face temperature is higher than T_Slag, then kill the particle and return

if(F_T(f,t)>T_Slag)
{
F_UDMI(f,t,4)+=flux/RHO_Slag;
C_UDMI(c0,t0,4)=F_UDMI(f,t,4);
return(PATH_ABORT);
}

//Add thickness to sintered layer if the face temperature is higher than T_Sinter, then kill particle and return

if(F_T(f,t)>T_Sinter)
{
F_UDMI(f,t,3)+=flux/RHO_Sinter;
C_UDMI(c0,t0,3)=F_UDMI(f,t,3);
return(PATH_ABORT);
}
//Add thickness to particulate layer if the surface temperature is lower than T_Sinter, then kill particle and return

F_UDMI(f,t,2)+=flux/RHO_Particulate;
C_UDMI(c0,t0,2)=F_UDMI(f,t,2);
return(PATH_ABORT);

}


DEFINE_DPM_INJECTION_INIT(Injection_Init,I)
{
thread_loop_f(t,d)
{

if(NNULLP(THREAD_STORAGE(t,SV_UDM_I)))
{

begin_f_loop(f,t);

{
L_Part=F_UDMI(f,t,2);
L_Sint=F_UDMI(f,t,3);
L_Slag=F_UDMI(f,t,4);
R_Part=L_Part/K_Particulate;
R_Sint=L_Sint/K_Sinter;
R_Slag=L_Slag/K_Slag;
R_Total=R_Part+R_Sint+R_Slag+R_BASE;
T_Surface=F_T(f,t);
cell_t c0;
Thread *t0;
c0=F_C0(f,t);
t0=F_C0_THREAD(f,t);
//Calculate interface temperatures between ash layers
T_PartSint_ifc=T_Cool+(T_Surface-T_Cool)*(R_Part+R_BASE)/R_Total;
T_SintSlag_ifc=T_Cool+(T_Surface-T_Cool)*(R_BASE+R_Part+R_Sint)/R_Total;
//Ckeck Melting and melt appropriate amount of sintered layer
if(T_SintSlag_ifc>T_Slag)
{

Y_Sint=R_Total*K_Sinter*(T_Slag-T_Cool)/(T_Surface-T_Cool)-K_Sinter*(R_Part+R_BASE);
Delta_L_Melt=L_Sint-Y_Sint;
F_UDMI(f,t,3)+=Y_Sint;
F_UDMI(f,t,4)+=L_Slag+Delta_L_Melt*RHO_Slag/RHO_Sinter;
C_UDMI(c0,t0,3)=F_UDMI(f,t,3);
C_UDMI(c0,t0,4)=F_UDMI(f,t,4);
}
//Check sintering and sinter appropriate amount of particulate layer
if(T_PartSint_ifc>T_Sinter);
{
Y_Part=R_Total*K_Particulate*(T_Sinter-T_Cool)/(T_Surface-T_Cool)-R_BASE*K_Particulate;
Delta_L_Sint=L_Part-Y_Part;
F_UDMI(f,t,2)+=Y_Part;
F_UDMI(f,t,3)+=L_Sint+Delta_L_Sint*RHO_Sinter/RHO_Particulate;
C_UDMI(c0,t0,2)=F_UDMI(f,t,2);
C_UDMI(c0,t0,3)=F_UDMI(f,t,3);
}
}
end_f_loop(f,t)
}
}

[pakk]

You change UDM3 in two places. Which one makes it negative?

The method for solving such problems is to make the code smaller, until you have the smallest code that gives the problem. Your code stoll sets UDM 0,1 and 2, which are irrelevant to the problem. Simplify it...

And your code is not complete, T_sinter and RHO_sinter are undefined.

[BW1211]

Originally Posted by pakk

You change UDM3 in two places. Which one makes it negative?

The method for solving such problems is to make the code smaller, until you have the smallest code that gives the problem. Your code stoll sets UDM 0,1 and 2, which are irrelevant to the problem. Simplify it...

And your code is not complete, T_sinter and RHO_sinter are undefined.

Thank you for your valuable reply. UDM3 is indeed used twice and this may be the problem. UDM3 in DEFINE_DPM_BC represents the thickness of the ash deposition, the value of each layer The thickness is not fixed. In DEFINE_DPM_INJECTION_INIT， UDM3 is also used as a variable to calculate the thickness of sinter layer. The particulate layer can transform into sinter layer and the sinter layer can also transform into the slag layer.
I think there was an error when the different thickness was converted to each other. T_sinter and RHO_Sinter are defined as global variables，UDM0,1,2 are necessary to characterize the deposition characteristics.
Thank you very much for your patience and valuable suggestions, I wish you a happy life.

[AlexanderZ]

in DEFINE_DPM_INJECTION_INIT you are using

L_Part=F_UDMI(f,t,2);
L_Sint=F_UDMI(f,t,3);
L_Slag=F_UDMI(f,t,4);

but I don't see where F_UDMI(f,t,2), .... are defined
so L_Part .... may have random values from the very beginning.

that could be the reason of your problem

[BW1211]

Originally Posted by AlexanderZ

in DEFINE_DPM_INJECTION_INIT you are using

L_Part=F_UDMI(f,t,2);
L_Sint=F_UDMI(f,t,3);
L_Slag=F_UDMI(f,t,4);

but I don't see where F_UDMI(f,t,2), .... are defined
so L_Part .... may have random values from the very beginning.

that could be the reason of your problem

Thank for your reply. Sorry for the unclear statement, but F_UDMI(f,t,2)...are defined in DEFINE_DPM_BC as below:
if(F_T(f,t)>T_Slag)
{
F_UDMI(f,t,4)+=flux/RHO_Slag;
C_UDMI(c0,t0,4)=F_UDMI(f,t,4);
return(PATH_ABORT);
}

//Add thickness to sintered layer if the face temperature is higher than T_Sinter, then kill particle and return

if(F_T(f,t)>T_Sinter)
{
F_UDMI(f,t,3)+=flux/RHO_Sinter;
C_UDMI(c0,t0,3)=F_UDMI(f,t,3);
return(PATH_ABORT);
}
//Add thickness to particulate layer if the surface temperature is lower than T_Sinter, then kill particle and return

F_UDMI(f,t,2)+=flux/RHO_Particulate;
C_UDMI(c0,t0,2)=F_UDMI(f,t,2);
return(PATH_ABORT);

[AlexanderZ]

how could you know that DEFINE_DPM_BC is executed before DEFINE_DPM_INJECTION_INIT

I'm not sure here, but anyway, for this expression

F_UDMI(f,t,4)+=flux/RHO_Slag;

you need previous value for summation, isn't it?
and previous value could be random

you may try to define explicitly values for

L_Part=F_UDMI(f,t,2);
L_Sint=F_UDMI(f,t,3);
L_Slag=F_UDMI(f,t,4);

for very first time step

[BW1211]

Originally Posted by AlexanderZ

how could you know that DEFINE_DPM_BC is executed before DEFINE_DPM_INJECTION_INIT

I'm not sure here, but anyway, for this expression

F_UDMI(f,t,4)+=flux/RHO_Slag;

you need previous value for summation, isn't it?
and previous value could be random

you may try to define explicitly values for

L_Part=F_UDMI(f,t,2);
L_Sint=F_UDMI(f,t,3);
L_Slag=F_UDMI(f,t,4);

for very first time step

Thank you for your valuable suggestion, I will have a try.

[pakk]

Originally Posted by BW1211

Thank you for your valuable reply. UDM3 is indeed used twice and this may be the problem.

Then investigate!

Remove one, do you still get negative results? Make the problem smaller, until it's just a few lines!

[BW1211]

Originally Posted by pakk

Then investigate!

Remove one, do you still get negative results? Make the problem smaller, until it's just a few lines!

Thank you very much. I will debug the udf code according to your suggestion.