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));
}
Show Posts