[jps115]

Hallo everyone,

I am a Master student in Transport Phenomena and for my thesis I am modeling the human respiratory system. In my thesis I am comparing lagrangian and Eulerian models for particle tracking.
I am trying to write a UDF that adds a Deposition flux for particles near the wall of my geometry in an Eulerian model.

I do this by first writing a DEFINE_ON_DEMAND udf that sets a UDMI to 0 for cells in my fluid and 1 to cells near the wall.

I then write a DEFINE_SOURCE udf that uses this UDMI value to see when to implement the flux. The UDF both compile and load with out error. But after the first iteration FLUENT crashes. I have atteched my code below.
The equation I am trying to implement is as follows:

Where is the convective velocity






If someone can help me out finding my mistake or can help me find a better approach to implement this flux it would be really appreciated.

#include "udf.h"
#include "flow.h"
#include "dpm.h"
#include "mem.h"
#include "sg.h"
#include "metric.h"

DEFINE_ON_DEMAND(selecting_the_cells)
{
	#if !RP_HOST
		real A[ND_ND];
		Domain *d;
		Thread *t, *tc, *t0;
		face_t f;
		cell_t c, c0;
		int Zone_ID = 13; /* Zone Id can be seen in the Boundary Conditions Panel, for me 13 = Wall */
		d=Get_Domain(1);
		tc=Lookup_Thread(d, Zone_ID); /* thread pointer of the wall */
		/* Loop through all the cell threads in the domain */
		thread_loop_c(t,d)
		{
			/* Loop through the cells in each cell thread */
			begin_c_loop(c,t)
			{
				C_UDMI(c,t,0)=0;
				C_UDMI(c,t,1)=0;
				C_UDMI(c,t,2)=0;
				C_UDMI(c,t,3)=0;
	
			}
			end_c_loop(c,t)
		}
		
		begin_f_loop(f,tc)
		{
			/* c0 and t0 identify the adjacent cell */
			c0=F_C0(f,tc);
			t0=THREAD_T0(tc);
			/* this loops over all cells adjacent to wall and lets the UDM = 1.0 */
			C_UDMI(c0,t0,0)=1.0;
			F_AREA(A,f,tc);
			C_UDMI(c0,t0,1)=A[0];
			C_UDMI(c0,t0,2)=A[1];
			C_UDMI(c0,t0,3)=A[2];
		}
		end_f_loop(f,tc)
	#endif
	Message("done with ON_DEMAND");
}

DEFINE_SOURCE(deposition_source, c, t, dS, eqn)
{
	real source;	
	real NV_VEC(v_c_vec);/* declaring vectors v_c_vec*/
	real NV_VEC(v_grad_v); /* declaring v*gradv as vector */
	real NV_VEC(diff); /* declaring diff as vector */
	real NV_VEC(psi_vec);/* declaring vectors psi */
	real NV_VEC(M_gravity); /* declaring vectors gravity*/
	real conv_flux, dif_flux;/* declaring value flux */
	real NV_VEC(A1), NV_VEC(A3); /* declaring vectors A */
	real grad_vel[3][3]; /* declaring array grad vel */
	real vel[3]; /* declaring array vel */
	float Um = 1.279279;/* Mean velocity at inlet m/s */
	float d1 = 0.006; /* Diameter of parent tube in m */
	float rho_p = 3413; /* density particle kg/m3 */
	double visc_f = 1.82e-5; /* kinematic viscosity */
	double d_p = 7.35e-6; /* particle diamter in m */
	double D = 1e-12; /* Mass Diffusivity m2/s (Stokes–Einstein equation) */
	NV_D(M_gravity, =, 0.0, 0.0, 9.81);
	real gravity = NV_MAG(M_gravity); /* M_gravity is the gravity vector M_gravity[0] = 0, M_gravity[1] = 0, M_gravity[2] = 9.81 */
	real St, Fr, Fr_inv, Pe, Pe_inv , A2[3], conc;
	St = rho_p*pow(d_p, 2)*Um/(18*visc_f*d1);
	Fr = pow(Um, 2)*gravity*d1;
	Pe = d1*Um/D;
	Pe_inv = pow(Pe, -1);
	Fr_inv = pow(Fr, -1);
	int i, j;
	conc = C_YI(c, t, 0) * rho_p;
	grad_vel[0][0] = C_DUDX(c, t), grad_vel[1][0] = C_DVDX(c, t), grad_vel[2][0] = C_DWDX(c, t);
	grad_vel[0][1] = C_DUDY(c, t), grad_vel[1][1] = C_DVDY(c, t), grad_vel[2][1] = C_DWDY(c, t);
	grad_vel[0][2] = C_DUDZ(c, t), grad_vel[1][2] = C_DVDZ(c, t), grad_vel[2][2] = C_DWDZ(c, t);

	vel[0] = C_U(c, t), vel[1] = C_V(c, t), vel[2] = C_W(c, t);
	NV_D(psi_vec, =, C_U(c,t),C_V(c,t),C_W(c,t));  /* defining psi in terms of velocity field */
	NV_V(A1, =, M_gravity);
	NV_S(A1, *=, Fr_inv);
	/* calculating vector matrix multiplication */
	for ( i = 0; i < 3; i++ )
	{
		A2[i] = 0;
		for ( j = 0; j < 3; j++ )
		A2[i] += grad_vel[i][j] * vel[j];
	}
	NV_D(A3, =, C_UDMI(c, t, 1), C_UDMI(c, t, 2), C_UDMI(c, t, 3));

	if (C_UDMI(c,t,0)>0.)
	{
	
		NV_D(v_grad_v, =, A2[0], A2[1], A2[2]); /* defining v_grad_v */
		NV_VS_VS(diff, =, A1, *, St, -, v_grad_v, *, St);
		NV_VV(v_c_vec, =, psi_vec, +, diff);
		NV_S(v_c_vec, *=, conc);
		conv_flux = NV_DOT(v_c_vec, A3);
		dif_flux = -Pe_inv*NV_DOT(C_YI_G(c,t,0), A3);
		C_UDMI(c, t, 4) += conv_flux + dif_flux;
		source = conv_flux + dif_flux;
		dS[eqn]= 0;
	}
	else
	{	
		C_UDMI(c, t, 4) = 0;
		source=  0;
		dS[eqn]=0.;
	}
	return source;
	
}

[jps115]

I noticed that "NV_DOT(C_YI_G(c,t,0), A3)" is the main problem in my code and cause fluent to crash. Reading the Fluent UDF manual it does state that C_YI_G(c,t,0) is a gradient vector. In my code A3 is defined as NV_D(A3, =, C_UDMI(c, t, 1), C_UDMI(c, t, 2), C_UDMI(c, t, 3)) with the C_UDMI being the vector components of the area vector. It seems as if the dot product is giving the problem. Sadly I am unable to understand why this is happening.

[AlexanderZ]

are you sure gradient is defined?
C_YI_G(c,t,0)

read ansys fluent customization manual -> Gradient (G) Vector Macros
manual says

you can
prevent the solver from freeing up memory by issuing the text command solve/set/expert and
then answering yes to the question, “Keep temporary solver memory from being freed?”

[jps115]

Hey Alexander,

Thank you for your response!
This does seem the case if I replace C_YI_G(c,t,0) with a defined vector it works fine. However solve/set/expert does not seem to be able to fix it. I will try a few things to get it to work. Probably if I run the simulation a few iteration before adding the UDF and setting this before solve/set/expert it will hopefully has some data for C_YI_G(c,t,0) and then run.

Thank you for your help. I will post an update if I am able to resolve it myself

[jps115]

I was able to identify the problem. I over looked that for the mass fraction gradient, it is written in the UDF manual that the mass fraction gradient is available only for the density-based solver and to use it in the pressure-based solver, I need to set the rpvar 'species/save-gradients? to #t. Hence, I typed in the text command the following

(rpsetvar 'species/save-gradients? #t)

This fixed everything. So in the end I should have read the manual more carefully. Thanks to everyone that helped me.