13

« **on:** April 25, 2012, 11:30:34 AM »
Dear piso

I am trying to do what u have mentioned in that thread. But my problem is when I try to write the location of each particle along with the ID of the tracked particle, it will give some erroneous result.

Bellow I am giving my UDF. I should know the exact location of each particle at any instant of time within the domain, so that with this information I can calculate the force.

Pl. reply.

#include "udf.h"

#define PI 3.14159 /*Assighned value for Pie*/

#define R_imp 1e-6 /*Radius of the Implant*/

#define mu_imp 2 /*Relative permeability of the Implant material*/

#define Tesla .6 /*Product of u0*H0*/

#define xfmp .4 /*weight fraction of froomagnetic material in MDCP*/

#define rhofmp 7850 /*Density of the ferromagnetic material*/

#define rhopol 950 /*Density of the binding polymer*/

#define r_mdcp 50e-9 /*radius of MDCP*/

#define Mfmp_s 1.735e6 /*Saturation Magnetization of Ferromagetic material*/

#define Bolt_const 1.3806503e-23 /*Boltzman constant*/

#define Temp 309 /*Absolute Temperature*/

real x,y,theta;

real Numr1,Dnmr1,Term1,Numr2,Dnmr2,Term2,Term3,Numr3,Numr4,Numr5,Numr6,Dnmr3,Dnmr4,Dnmr5,Dnmr6;

real Term_sq,C1,Del_Phi_del_x,Del_Phi_del_y,Del2_Phi_del_x2,Del2_Phi_del_y2,Del2_Phi_del_xy;

real B_x,B_y,mod_B;

real Del_Bx_delx,Del_Bx_dely,Del_By_delx,Del_By_dely;

real wfmp_denom,wfmp_numr,wfmp;

real V_p,beta_numr,beta_dnmr,beta1,beta,Lan_beta;

real moment1,moment2,moment3,moment4,moment_x,moment_y;

real Fmx1,Fmx2,Fmx,Fmy1,Fmy2,Fmy,bforce;

float x1,y1,z1;

DEFINE_DPM_BODY_FORCE(particle_body_force,p,i)

{

FILE *pf;

theta=0;

x=P_POS(p)[0];

y=P_POS(p)[1];

Term_sq=(x*x+y*y);

x1=p->state.pos[0];

y1=p->state.pos[1];

z1=p->state.pos[2];

if(NULL == (pf = fopen("particle_position.txt","a")))

Error("Could not open file for append!\n");

fprintf(pf,"%d\t%e\t%e\n",p,x1,y1);

fclose(pf);

/*calculation of Term1 in expression of Phi*/

Numr1=mu_imp-1;

Dnmr1=mu_imp+1;

Term1=Numr1/Dnmr1;

/*calculation of Term2 in expression of Phi*/

Term2=R_imp*R_imp;

C1=Term1*Term2;

/*calculation of Dphi/dx*/

Numr2=((cos(theta)*(y*y-x*x))-(2*x*y*sin(theta)));

Dnmr2=Term_sq*Term_sq;

Del_Phi_del_x=C1*Numr2/Dnmr2;

/*calculation of Dphi/dy*/

Numr3=((sin(theta)*(x*x-y*y))-(2*x*y*cos(theta)));

Dnmr3=Term_sq*Term_sq;

Del_Phi_del_y=C1*Numr3/Dnmr3;

/*calculation of D2phi/dx2*/

Numr4=((x*cos(theta)*(x*x-3*y*y))+(y*sin(theta)*(3*x*x-y*y)));

Dnmr4=Term_sq*Term_sq*Term_sq;

Del2_Phi_del_x2=2*Tesla*C1*Numr4/Dnmr4;

/*calculation of D2phi/dy2*/

Numr5=((x*cos(theta)*(-x*x+3*y*y))+(y*sin(theta)*(-3*x*x+y*y)));

Dnmr5=Term_sq*Term_sq*Term_sq;

Del2_Phi_del_y2=2*Tesla*C1*Numr5/Dnmr5;

/*calculation of D2phi/dxdy*/

Numr6=((y*cos(theta)*(3*x*x-y*y))+(x*sin(theta)*(3*y*y-x*x)));

Dnmr6=Term_sq*Term_sq*Term_sq;

Del2_Phi_del_xy=2*Tesla*C1*Numr6/Dnmr6;

/*calculation of B*/

B_x=(Tesla)*cos(theta)-(Tesla)*Del_Phi_del_x;

B_y=(Tesla)*sin(theta)-(Tesla)*Del_Phi_del_y;;

mod_B=sqrt((B_x*B_x)+(B_y*B_y));

/*calculation of volume fraction*/

wfmp_denom=(xfmp+((1-xfmp)*(rhofmp/rhopol)));

wfmp_numr=xfmp;

wfmp=wfmp_numr/wfmp_denom;

/*calculation of volume of MDCP*/

V_p=(4/3)*PI*(r_mdcp*r_mdcp*r_mdcp);

/*calculation of Beta*/

beta_numr=wfmp*V_p*Mfmp_s;

beta_dnmr=Bolt_const*Temp;

beta1=beta_numr/beta_dnmr;

beta=beta1*mod_B;

/*calculation of langevin function*/

Lan_beta=1/(tanh(beta))-1/beta;

/*calculation of magnetic moment*/

moment1=wfmp*V_p*Mfmp_s;

moment2=Lan_beta;

moment3=mod_B;

moment4=(moment1*moment2)/moment3;

moment_x=moment4*B_x;

moment_y=moment4*B_y;

if(i==0.)

{

Fmx1=moment_x*Del2_Phi_del_x2;

Fmx2=moment_y*Del2_Phi_del_xy;

Fmx=Fmx1+Fmx2;

bforce=-Fmx;

}

else if(i==1.)

{

Fmy1=moment_x*Del2_Phi_del_xy;

Fmy2=moment_y*Del2_Phi_del_y2;

Fmy=Fmy1+Fmy2;

bforce=-Fmy;

}

return(bforce/P_MASS(p));

}