
[Sponsors] 
Deposition flux for particles in an Eulerian model 

LinkBack  Thread Tools  Search this Thread  Display Modes 
February 17, 2022, 04:11 
Deposition flux for particles in an Eulerian model

#1 
New Member
Jeroen de Graaff
Join Date: Feb 2022
Posts: 9
Rep Power: 4 
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. Code:
#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.82e5; /* kinematic viscosity */ double d_p = 7.35e6; /* particle diamter in m */ double D = 1e12; /* 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; } 

February 18, 2022, 06:53 

#2 
New Member
Jeroen de Graaff
Join Date: Feb 2022
Posts: 9
Rep Power: 4 
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.


February 20, 2022, 21:37 

#3  
Senior Member
Alexander
Join Date: Apr 2013
Posts: 2,363
Rep Power: 34 
are you sure gradient is defined?
C_YI_G(c,t,0) read ansys fluent customization manual > Gradient (G) Vector Macros manual says Quote:
__________________
best regards ****************************** press LIKE if this message was helpful 

February 21, 2022, 05:15 

#4 
New Member
Jeroen de Graaff
Join Date: Feb 2022
Posts: 9
Rep Power: 4 
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 

February 21, 2022, 06:08 

#5 
New Member
Jeroen de Graaff
Join Date: Feb 2022
Posts: 9
Rep Power: 4 
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 densitybased solver and to use it in the pressurebased solver, I need to set the rpvar 'species/savegradients? to #t. Hence, I typed in the text command the following
(rpsetvar 'species/savegradients? #t) This fixed everything. So in the end I should have read the manual more carefully. Thanks to everyone that helped me. 

Tags 
deposition, eulerian, flux calculation, udf 
Thread Tools  Search this Thread 
Display Modes  


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  STARCCM+  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 