CFD Online Discussion Forums - Deposition flux for particles in an Eulerian model

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.

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.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;

        

}