CFD Online Discussion Forums

I wrote the udf after going through Fluent user guide. Can any one go through it and let me know if I am right in calculating force and moment vectors ..

/*Fluent Six-Degree of freedom User Defined Function (6-DOF udf).*/

#include "udf.h"

/*Define mass */

#DEFINE mass 50

/*Define Moment of Inertia about x-axis*/

#DEFINE Ix 5

/*Define Moment of Inertia about y-axis*/

#DEFINE Iy 5

/*Define Moment of Inertia about z-axis*/

#DEFINE Iz 5

/*Define Center of Gravity*/

#DEFINE cg_x 0

#DEFINE cg_y 0

#DEFINE cg_z 0

static real vx = 0.0;

static real vy = 0.0;

static real vx = 0.0;

static real wx = 0.0;

static real wy = 0.0;

static real wz = 0.0;

DEFINE_CG_MOTION(block,dt,vel,omega,time,dtime)

{ Thread *t;

face_t f;

Node *node;

int n;

real NV_VEC(A);

real Fx, Fy, Fz, Mx, My, Mz;

real dvx, dvy, dvz, dwx, dwy, dwz;

/* reset velocities */

NV_S(vel, =, 0.0);

NV_S(omega, =, 0.0);

if (!Data_Valid_P())

return;

/* get the thread pointer for which this motion is defined */

t = DT_THREAD(dt);

/* compute pressure force on body by looping through all faces */

Fx = 0.0;

Fy = 0.0;

Fx = 0.0;

Mx = 0.0;

My = 0.0;

Mz = 0.0;

begin_f_loop(f,t)

{

f_node_loop(f,t,n)

{

F_AREA(A,f,t);

/* x,y,z - components of Force

x-component of force = Pressure on the face * x-component of Area vector

*/

Fx += F_P(f,t)*A[0];

Fy += F_P(f,t)*A[1];

Fz += F_P(f,t)*A[2];

/* Moment (M) = vector product of Force and Displacement vector (F X r);

F = Fx i + Fy j +Fz k

r = (x - cg_x)i + (y- cg_y)j + (z - cg_z)k

*/

Mx += F_P(f,t)*(A[1]*(NODE_Z(node) - cg_z)-A[2]*(NODE_Y(node) - cg_y));

My += F_P(f,t)*(A[2]*(NODE_X(node) - cg_x)-A[0]*(NODE_Z(node) - cg_z));

Mx += F_P(f,t)*(A[0]*(NODE_Y(node) - cg_y)-A[1]*(NODE_X(node) - cg_x));

}

} end_f_loop(f,t)

/* compute change in velocity, i.e., dv = F * dt / mass

velocity update using explicit Euler formula */

dvx = dtime * Fx / mass;

vx +=dvx;

dvy = dtime * Fy / mass;

vy +=dvy;

dvz = dtime * Fz / mass;

vz +=dvz;

dwx = dtime * Mx / Ix;

wx +=dwx;

dwy = dtime * My / Iy;

wy +=dwy;

dwz = dtime * Mz / Iz;

wz +=dwz;

Message ("time = %f\n", time);

Message ("x_force = %f, y_force = %f, z_force = %f\n", Fx, Fy, Fz); Message ("x_vel = %f, y_vel = %f, z_vel = %f\n", vx, vy, vz);

Message ("x_moment = %f, y_moment = %f, z_moment = %f\n", Mx, My, Mz); Message ("x_omega = %f, y_omega = %f, z_omega = %f\n", wx, wy, wz);

/* set x,y,z-component of velocity */

vel[0] = vx;

vel[1] = vy;

vel[2] = vz;

/* set x,y,z-component of angular velocity */

omega[0]= wx;

omega[1]= wy;

omega[2]= wz;

}