CFD Online Discussion Forums - udf for valve closing a pipe using dynamic mesh

hi everyone. I want to model a valve closing in the start of a pipe. there is a tutorial in fluent help, describing this problem, but there is a difference between this tutorial and my problem. in the tutorial the valve moves horizontally to close the pipe but I want the vale to move vertically to close the pipe. I mean it should enter the pipe from top or below of the pipe. my question is that how should I change the udf to make sensible for my problem. the following is the picture of problem description and udf used for it in the fluent tutorial help. I really be thankful if anyone could help me change it for my own problem.

# include "udf.h"
#
include "dynamesh_tools.h"

/********************** UDF Explanation Starts ***************************
This udf is used to solve transient valve solidy body Fluid Structure
Interaction problem. A more sophisticated version for more general 1dof
problem is also availabe upon request.

Below are some explanations about how to use the udf:

valveid[NO_OF_VALVES][NO_OF_ZONES]={{22, 23, -1}}: valve face
zone list, in this example, the valve consists of face zone 21 and 23,
-1 is a flag and so keep it. These face zones will be used to calculate
the force on the valve.

lift_min[NO_OF_VALVES]={0.0005}: minimum lift.
lift_max[NO_OF_VALVES]={0.006}: maximum lift
rest_conts[NO_OF_VALVES]={0.5}: restitution factor
mass[NO_OF_VALVES]={0.02}: valve mass
stiffness[NO_OF_VALVES]={2000}: valve spring stiffness
stretch_at_closed[NO_OF_VALVES]={0.002}: the valve stretching
length when the valve is closed. The spring force is supposed to pull
the valve from the max lift location to min lift location.

current_vel_mag[NO_OF_VALVES]: Initial velocity magnitude. The value is negative
if the valve moves toward positive x direction (yes,
it is a little strange).
INITIAL_LIFT: Initial valve lift

How to use the udf:

- Modify the user input part of the udf.
- Build the lib and hook the reader/writer udf
- Use the udf for the valve motion

Limitations:

- The valve can consist of up to 10 face zones.
- Single valve only.
- The valve has to move in the x direction with min lift location to the left
of the max lift location.

Written by : Xiao Hu (Fluent Inc)
Last updated : 02/24/2006

********************** UDF Explanation Ends *****************************/

# define DEBUG 1
# define NO_OF_VALVES 1
# define NO_OF_ZONES 11

/*************************** User Input Starts *****************************/

# define INITIAL_LIFT 10e-3

static int valveid[NO_OF_VALVES][NO_OF_ZONES]={{7, -1}};

static real lift_min[NO_OF_VALVES]={1e-3};
static real lift_max[NO_OF_VALVES]={1000e-3};
static real mass[NO_OF_VALVES]={0.2};
static real stiffness[NO_OF_VALVES]={1000};
static real stretch_at_closed[NO_OF_VALVES]={11e-3};

static real rest_conts[NO_OF_VALVES]={0};
static real current_vel_mag[NO_OF_VALVES]={0};

/*************************** User Input Ends *******************************/

static real axis[NO_OF_VALVES][ND_ND]={{-1, 0}}; /* normalized */
static real gravity_direction[ND_ND]={0, 1}; /* normalized */
static real r_rp_closed[NO_OF_VALVES][ND_ND]={{0,0}};
static real cur_r_rp[NO_OF_VALVES][ND_ND]={{INITIAL_LIFT, 0}};
static real previous_time[NO_OF_VALVES]={0};

static void f_valve(int valveNo, void *dt, real *cg_vel, real *cg_omega, real time, real dtime)
{
#if !RP_HOST
real tmp[ND_ND], dv, current_vel[ND_ND], CG[ND_ND], force[3], moment[3], stretch;
real aero_force[ND_ND], aero_force_axis, spring_force, gravity_force, net_force,
r_rp_new[NO_OF_VALVES][ND_ND];
int i;
Thread * tf;
Domain * domain;

/************************************************** ****************/

static real cg_vel_saved[NO_OF_VALVES][ND_ND];

/************************************************** ****************/

/* Do the calculation if the new time step */

if(fabs(previous_time[valveNo]-time)>0.2*dtime)

{
/* reset velocities */

NV_S (cg_vel, =, 0.0);

NV_S (cg_omega, =, 0.0);

/* Check to see if there is data */

if (!Data_Valid_P ())

{

Message0("\n\nNo data->No mesh motion!!!\n\n");

return;

}

/*Calculate force*/

domain = THREAD_DOMAIN (DT_THREAD ((Dynamic_Thread *)dt));

i=0;

NV_S(aero_force,=,0);

while(valveid[valveNo][i]>=0)

{

tf=Lookup_Thread(domain, valveid[valveNo][i]);

NV_S (CG, =, 0.0);

/* This function needs size of 3 for force and moment. It works in parallel and
the common arguments are the of the same type as your arguments.
In addition you have to pass the domain as first arg and a boolean as
last argument. This boolean has to be TRUE if you also call the function
on the host else it has to be FALSE. */

Compute_Force_And_Moment (domain, tf, CG, force, moment, FALSE);

NV_V(aero_force,+=,force);

i++;

}

aero_force_axis=NV_DOT(aero_force, axis[valveNo]);

NV_VV(tmp,=,r_rp_closed[valveNo],-,cur_r_rp[valveNo]);

stretch = (stretch_at_closed[valveNo]+NV_DOT(tmp,axis[valveNo]));

spring_force=stiffness[valveNo]*stretch;

gravity_force = 9.81*mass[valveNo]*NV_DOT(gravity_direction, axis[valveNo]);

net_force=spring_force+aero_force_axis+gravity_for ce;

dv=net_force/mass[valveNo]*dtime;

/* Calculate the C.G location and velocity if it does not hit the boundary */

NV_VS(current_vel,=,axis[valveNo],*,current_vel_mag[valveNo]);

NV_VS(tmp,=,current_vel,*,dtime);

NV_VV(r_rp_new[valveNo],=,cur_r_rp[valveNo],+,tmp);

/* Update velocity */

current_vel_mag[valveNo]+=dv;

/* debug info */

#if DEBUG

Message0("\n\n*********************** VALVE INFO ***************************\n");

Message0("\nCurrent valve lift =%10.3e\n", -NV_DOT(cur_r_rp[valveNo], axis[valveNo]));

Message0("\nCurrent valve velocity (negative if pointing to positive x) =%10.3e\n", current_vel_mag[valveNo]-dv);

Message0("\naero force =%10.3e\n", aero_force[0]);

Message0("\n(stretching at closed, stretching, force)=(%10.3e, %10.3e, %10.3e)\n",

stretch_at_closed[valveNo], stretch, -spring_force);

Message0("\n(net_force, spring force, aero force)=(%10.3e, %10.3e, %10.3e)\n",
-net_force, -spring_force, -aero_force_axis);

#endif

/* if it hits the lift_min boundary then it stays at lift_min*/

NV_VV(tmp,=,r_rp_closed[valveNo],-,r_rp_new[valveNo]);

if(NV_DOT(tmp,axis[valveNo])<1.0000001*(lift_min[valveNo]))

{

NV_V_VS(r_rp_new[valveNo],=,r_rp_closed[valveNo],-,axis[valveNo],*,lift_min[valveNo]);

current_vel_mag[valveNo]=-rest_conts[valveNo]*fabs(current_vel_mag[valveNo]);

#if DEBUG

Message0("\n Valve hits the min valve lift\n");

#endif

}

/* if it hits the lift_max boundary then it stays at lift_max*/

NV_VV(tmp,=,r_rp_closed[valveNo],-,r_rp_new[valveNo]);

if(NV_DOT(tmp,axis[valveNo])>0.99999999*(lift_max[valveNo]))

{

NV_V_VS(r_rp_new[valveNo],=,r_rp_closed[valveNo],-,axis[valveNo],*,lift_max[valveNo]);

current_vel_mag[valveNo]=rest_conts[valveNo]*fabs(current_vel_mag[valveNo]);

#if DEBUG

Message0("\n Valve hits the max valve lift\n");

#endif

}

/* set valve velocity */

NV_VV(tmp,=,r_rp_new[valveNo],-,cur_r_rp[valveNo]);

NV_VS(cg_vel,=,tmp,/,dtime);

/* Update location and velocity */

NV_V(cur_r_rp[valveNo],=,r_rp_new[valveNo]);

NV_V(cg_vel_saved[valveNo],=,cg_vel);

previous_time[valveNo]=time;

/* debug info */

#if DEBUG

Message0("\nNext valve velocity (negative if pointing to positive x)=%11.3e\n", current_vel_mag[valveNo]);

Message0("\nNext valve lift =%10.3e\n", -NV_DOT(r_rp_new[valveNo], axis[valveNo]));

Message0("\n*********************** VALVE INFO ***************************\n\n");

#endif

}

else

{

NV_V(cg_vel,=,cg_vel_saved[valveNo]);

}

#endif

node_to_host_real(current_vel_mag, NO_OF_VALVES);

node_to_host_real(cur_r_rp[0], NO_OF_VALVES*ND_ND);

node_to_host_real(previous_time, NO_OF_VALVES);

}

DEFINE_CG_MOTION(valve, dt, cg_vel, cg_omega, time, dtime)

{

f_valve(0, dt, cg_vel, cg_omega, time, dtime);

node_to_host_real(cg_vel,ND_ND);

node_to_host_real(cg_omega,ND_ND);

}

static void write_data(FILE *fp)

{

int i, j;

for(i=0; i<NO_OF_VALVES; i++)

{

fprintf(fp, "%e ", current_vel_mag[i]);

}

fprintf(fp, "\n");

for(i=0; i<NO_OF_VALVES; i++)

{

for(j=0; j<ND_ND; j++)

{

fprintf(fp, "%e ", cur_r_rp[i][j]);

}

fprintf(fp, "\n");

}

fprintf(fp, "\n");

for(i=0; i<NO_OF_VALVES; i++)

{

fprintf(fp, "%e ", previous_time[i]);

}

}

static void read_data(FILE * fp)

{

int i, j;

for(i=0; i<NO_OF_VALVES; i++)

{

#if RP_DOUBLE

fscanf(fp, "%le", current_vel_mag+i);

#else

fscanf(fp, "%e", current_vel_mag+i);

#endif

}

for(i=0; i<NO_OF_VALVES; i++)

{

for(j=0; j<ND_ND; j++)

{

#if RP_DOUBLE

fscanf(fp, "%le", cur_r_rp[i]+j);

#else

fscanf(fp, "%e", cur_r_rp[i]+j);

#endif

}

}

for(i=0; i<NO_OF_VALVES; i++)

{

#if RP_DOUBLE

fscanf(fp, "%le", previous_time+i);

#else

fscanf(fp, "%e", previous_time+i);

#endif

}

}

DEFINE_RW_FILE(writer, fp)

{

Message0("Writing UDF data to data file...\n");

#if PARALLEL

#if RP_HOST

write_data(fp);

#endif
#else

write_data(fp);

#endif

}

DEFINE_RW_FILE(reader, fp)

{

Message0("Reading UDF data from data file...\n");

#if PARALLEL

#if RP_HOST

read_data(fp);

#endif

#else

read_data(fp);

#endif

host_to_node_real(current_vel_mag, NO_OF_VALVES);

host_to_node_real(cur_r_rp[0], NO_OF_VALVES*ND_ND);

host_to_node_real(previous_time, NO_OF_VALVES);

}