CFD Online Logo CFD Online URL
Home > Forums > Software User Forums > ANSYS > FLUENT > Fluent UDF and Scheme Programming

Deposition flux for particles in an Eulerian model

Register Blogs Members List Search Today's Posts Mark Forums Read

LinkBack Thread Tools Search this Thread Display Modes
Old   February 17, 2022, 04:11
Default Deposition flux for particles in an Eulerian model
New Member
Jeroen de Graaff
Join Date: Feb 2022
Posts: 9
Rep Power: 4
jps115 is on a distinguished road
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:
\vec{vc} = \vec{v} + St(Fr^{-1}\vec{g}-\vec{v}\bullet\nabla\vec{v})
Where \vec{vc} is the convective velocity

F_{dep} = F^c|_{w}+F^d|_{w}

F^c|_{w} = \int c \vec{vc}\bullet\vec{d}S|_{w}
F^d|_{w} = -Pe^{-1} \int \nabla c \bullet\vec{d}S|_{w}

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"

	#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 */
		tc=Lookup_Thread(d, Zone_ID); /* thread pointer of the wall */
		/* Loop through all the cell threads in the domain */
			/* Loop through the cells in each cell thread */
			/* c0 and t0 identify the adjacent cell */
			/* this loops over all cells adjacent to wall and lets the UDM = 1.0 */
	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;
		C_UDMI(c, t, 4) = 0;
		source=  0;
	return source;
jps115 is offline   Reply With Quote

Old   February 18, 2022, 06:53
New Member
Jeroen de Graaff
Join Date: Feb 2022
Posts: 9
Rep Power: 4
jps115 is on a distinguished road
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.
jps115 is offline   Reply With Quote

Old   February 20, 2022, 21:37
Senior Member
Join Date: Apr 2013
Posts: 2,363
Rep Power: 34
AlexanderZ will become famous soon enoughAlexanderZ will become famous soon enough
are you sure gradient is defined?

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?”
best regards

press LIKE if this message was helpful
AlexanderZ is offline   Reply With Quote

Old   February 21, 2022, 05:15
New Member
Jeroen de Graaff
Join Date: Feb 2022
Posts: 9
Rep Power: 4
jps115 is on a distinguished road
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 is offline   Reply With Quote

Old   February 21, 2022, 06:08
New Member
Jeroen de Graaff
Join Date: Feb 2022
Posts: 9
Rep Power: 4
jps115 is on a distinguished road
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.
jps115 is offline   Reply With Quote


deposition, eulerian, flux calculation, udf

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On

Similar Threads
Thread Thread Starter Forum Replies Last Post
Wind turbine simulation Saturn CFX 58 July 3, 2020 01:13
Eulerian model theory hamed.majeed FLUENT 7 September 4, 2018 07:16
How to simulate the eulerian multiphase model about particle jhlee9622 STAR-CCM+ 2 November 24, 2016 11:37
Superlinear speedup in OpenFOAM 13 msrinath80 OpenFOAM Running, Solving & CFD 18 March 3, 2015 05:36
An error has occurred in cfx5solve: volo87 CFX 5 June 14, 2013 17:44

All times are GMT -4. The time now is 01:01.