把udf贴给大家吧!用法自己找书看看!- /**********************************************************************************************************
- UDF for vane pump application
- Mesh requirements:
- - Rotation must be long Z axis.
- - The origin must be placed on the cylindrical axis of rotating shaft or the SMALLER circle.
- - The large circle (if circle is used) must be placed to the left of the small circle
- and its center must have zero y coordinate.
- - Initially, one gap must be on the positive x location.
- (The first three are just orientation. The last one is required because the udf assumes this initial
- position for the angle alfa calculation.)
- - All the pump core needs to be a single cell zone and all gaps need to be a single
- cell zone.
- How to use the udf:
- - Modify the input file (input.txt) and place in working directory.
- The input file description is as follows;
- 0 - Outer hape type (0=circle, 1=user-defined)
- 0 - Inner shape type (0=circle, 1=user-defined, 2=special)
- 8 - Number of vanes
- 6 11 - (Core ID and Gap ID)
- 25e-3 1.75e-3 0.1e-3 (inner circle radius, half vane width and gap size)
- 27.5e-3 2e-3 (outer radius and offset - if it is a circle)
- - If outer profile is not a circle, then this profile must be created and stored in a file
- called data_outer.txt. This file must also be placed in the working directory and has the
- following format;
- 320 - Number of points
- 0.0275 0 - x and y coordinates of each point required to create the profile
- 0.027497851 0.0005
- 0.027491405 0.001
- 0.027480657 0.0015
- 0.027465603 0.002
- - Build the UDF library.
- - Load the mesh (The mesh has to satisfy the requirements above).
- - Turn on Dynamic Mesh and In-Cylinder.
- - Specify RPM under In-Cylinder tab. Under the same tab, use 360 for Crank period and 0.25 for
- Crank Angle Step Size (this determines time step size).
- - Hook the UDF with fluid-pump and fluid-gap.
- - Hook the initialization udf.
- - Hook the reader and writer udf.
- - Define one UDM.
- - Define grid interfaces.
- - Define other parameters for transient flow (turbulence, PISO etc.).
- - Initialize the flow (you will need to intialize the flow even for just mesh motion)
- Note:
- - Need to read both cas/dat files even for just mesh motion
- - Cannot initialize the flow in the middle of a run and then continue.
- Written by: Xiao Hu ([url=mailto:xh@fluent.com]xh@fluent.com[/url]), and Laz Foley ([url=mailto:lnf@fluent.com]lnf@fluent.com[/url])
- Last Updated: 1/20/2004
- **********************************************************************************************************/
- # include "udf.h"
- # include "dynamesh_tools.h"
- # define RPM RP_Get_Real("dynamesh/in-cyn/crank-rpm")
- # define noDEBUG
- /************************************* User inputs start **********************************************/
- /* extern real special_inner_radius(real z); */
- /* Number of vanes, core ID and gap ID */
- int N_VANE, CORE_ID, GAP_ID;
- /* Outer circle radius, inner circle radius, half vane width, gap size and offset */
- real R, r, HALF_VANE_WIDTH, GAP, DELTA;
- /* Outer contour shape */
- enum define_shape
- {
- circle, user_defined, special
- }OUTER_SHAPE, INNER_SHAPE;
- /************************************* User inputs end ************************************************/
- enum shrink
- {
- w_shrink, wo_shrink
- };
- struct shape
- {
- enum define_shape shape_type;
- enum shrink shrink_y_or_n;
- real center[2];
- real radius;
- real (*data)[2];
- int number_of_data;
- }pump_core_outer_shape, pump_core_inner_shape, pump_gap_outer_shape, pump_gap_inner_shape;
- enum UDM
- {
- UDM_sector
- };
- static real my_atan2(real y, real x)
- {
- real result;
- result=atan2(y,x);
- if(atan2(y,x)<0)
- result = atan2(y,x) + 2*M_PI;
- return result;
- }
- /* Find the intersection t1 and t2 for two lines (see notes pg 5)*/
- void find_t(real xa, real ya, real xb, real yb,
- real xc, real yc, real xd, real yd,
- real *t1, real *t2, int * flag)
- {
- real k;
- if((xd-xc)!=0)
- {
- k = (yd-yc)/(xd-xc);
- if((yb-ya)-k*(xb-xa)==0)
- {
- (*flag) = 0;
- return;
- }
- (*flag) = 1;
- (*t1) = (-k*(xc-xa)+(yc-ya))/((yb-ya)-k*(xb-xa));
- (*t2) = ((xb-xa)*(*t1)-(xc-xa))/(xd-xc);
- }
- else
- {
- if((xb-xa)==0)
- {
- (*flag) = 0;
- return;
- }
- if(yd==yc)
- {
- Message0("\nc and d are the same points-aborting!!!\n");
- exit(0);
- }
- (*flag) = 2;
- (*t1) = (xc-xa)/(xb-xa);
- (*t2) = ((yb-ya)*(*t1)-(yc-ya))/(yd-yc);
- }
- return;
- }
- /* Find the intersection (see notes pg5) */
- void find_intersection(real (*data)[2], int number_of_data, real xc, real yc, real xd, real yd,
- real *x, real *y)
- {
- int i, flag;
- real pt1[2], pt2[2], t1, t2;
- for(i=0; i<number_of_data; i++)
- {
- if(i<number_of_data-1)
- {
- pt1[0] = data[i][0];
- pt1[1] = data[i][1];
- pt2[0] = data[i+1][0];
- pt2[1] = data[i+1][1];
- }
- else
- {
- pt1[0] = data[i][0];
- pt1[1] = data[i][1];
- pt2[0] = data[0][0];
- pt2[1] = data[0][1];
- }
- find_t(pt1[0], pt1[1], pt2[0], pt2[1], xc, yc, xd, yd, &t1, &t2, &flag);
- /*Message0("\n%7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f flag=%3d\n",
- pt1[0], pt1[1], pt2[0], pt2[1], xc, yc, xd, yd, t1, t2, flag); */
- if(((flag==1)||(flag==2))&&(t1<=1)&&(t1>=0)&&(t2>=0))
- {
- (*x) = (1-t1)*pt1[0] + t1*pt2[0];
- (*y) = (1-t1)*pt1[1] + t1*pt2[1];
- /* Message0("\n%7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f \n", pt1[0], pt1[1], pt2[0], pt2[1], xc, yc, xd, yd);
- Message0("\nt1=%5.2f t2=%5.2f x=%5.2f y=%5.2f\n", t1, t2, *x, *y); */
- return;
- }
- }
- }
- /* Function used to calcuate op for the gap. Please refer to note page 4 for the definition. */
- static void f_Op_gap(real * xop, real * yop, real xa, real ya, real time, real sector)
- {
- real alfa, h;
- alfa = time*RPM*M_PI/30+sector*2*M_PI/N_VANE;
- h=cos(alfa)*ya-sin(alfa)*xa;
- if(fabs(cos(alfa))>0.1)
- {
- *xop = 0.5*xa;
- *yop = (sin(alfa)*0.5*xa+h)/cos(alfa);
- }
- else
- {
- *xop = (cos(alfa)*0.5*ya-h)/sin(alfa);
- *yop = 0.5*ya;
- }
- return;
- }
- /* Function used to calcuate op. Please refer to note page 3 for the definition. */
- static void f_Op_core(real * x, real * y, real time, real sector)
- {
- real alfa;
- alfa = time*RPM*M_PI/30+sector*2*M_PI/N_VANE;
- *x = HALF_VANE_WIDTH*(cos(alfa)+cos(alfa+2*M_PI/N_VANE))/sin(2*M_PI/N_VANE);
- *y = HALF_VANE_WIDTH*(sin(alfa)+sin(alfa+2*M_PI/N_VANE))/sin(2*M_PI/N_VANE);
- return;
- }
- /* Generic function to calcualte the distance between point a and c. Please refer to note page 2.
- The inputs are points a, b, center of the circle (x0, y0), and the circle radius. The h and H
- for this application are all from this function.
- */
- static real f_h(real xa, real ya, real xb, real yb, real z_coord, struct shape * contour_shape)
- {
- if (contour_shape->shape_type==circle)
- {
- real t1, t2, a, b, c, t, x0, y0, radius;
- x0 = contour_shape->center[0];
- y0 = contour_shape->center[1];
- radius = contour_shape->radius;
- a=pow(xb-xa,2)+pow(yb-ya,2);
- b=2*(xa-x0)*(xb-xa)+2*(ya-y0)*(yb-ya);
- c=pow(xa-x0,2)+pow(ya-y0,2)-radius*radius;
- t1 = (-b+sqrt(b*b-4*a*c))/(2*a);
- t2 = (-b-sqrt(b*b-4*a*c))/(2*a);
- if (t1>0&&t2<0)
- t = t1;
- else if (t2>0&&t1<0)
- t = t2;
- else
- {
- Message("\nSomething wrong with f_h. t1=%7.2f t2=%7.2f", t1, t2);
- Message("\nxa=%11.3e ya=%11.3e xb=%11.3e yb=%11.3e x0=%10.3e y0=%10.3e radius=%10.3e", xa, ya, xb, yb, x0, y0, radius);
- t = 1;
- }
- return t*sqrt((pow(yb-ya,2)+pow(xb-xa,2)));
- }
- else if (contour_shape->shape_type==user_defined)
- {
- real x, y;
- find_intersection(contour_shape->data, contour_shape->number_of_data, xa, ya, xb, yb, &x, &y);
- if(contour_shape->shrink_y_or_n==w_shrink)
- {
- return sqrt(pow(y-ya,2)+pow(x-xa,2))-GAP;
- }
- else
- {
- return sqrt(pow(y-ya,2)+pow(x-xa,2));
- }
- }
- else if (contour_shape->shape_type==special)
- {
- real t1, t2, a, b, c, t, x0, y0, radius;
- x0 = 0;
- y0 = 0;
- radius = 1;
- /* radius = special_inner_radius(z_coord); */
- a=pow(xb-xa,2)+pow(yb-ya,2);
- b=2*(xa-x0)*(xb-xa)+2*(ya-y0)*(yb-ya);
- c=pow(xa-x0,2)+pow(ya-y0,2)-radius*radius;
- t1 = (-b+sqrt(b*b-4*a*c))/(2*a);
- t2 = (-b-sqrt(b*b-4*a*c))/(2*a);
- if (t1>0&&t2<0)
- t = t1;
- else if (t2>0&&t1<0)
- t = t2;
- else
- {
- Message("\nSomething wrong with f_h. t1=%7.2f t2=%7.2f", t1, t2);
- Message("\nxa=%11.3e ya=%11.3e xb=%11.3e yb=%11.3e x0=%10.3e y0=%10.3e radius=%10.3e", xa, ya, xb, yb, x0, y0, radius);
- t = 1;
- }
- return t*sqrt((pow(yb-ya,2)+pow(xb-xa,2)));
- }
- else
- {
- Message0("\nwrong type-aborting!!!\n");
- exit(0);
- }
- }
- static real f_h_core(real * op, real *ap, real time)
- {
- return f_h(op[0], op[1], ap[0], ap[1], ap[2], &pump_core_inner_shape);
- }
- static real f_H_core(real * op, real *ap, real time)
- {
- real not_used_var=0;
- return f_h(op[0], op[1], ap[0], ap[1], not_used_var, &pump_core_outer_shape);
- }
- static real f_h_gap(real * op, real *ap, real time)
- {
- real not_used_var=0;
- return f_h(op[0], op[1], ap[0], ap[1], not_used_var, &pump_gap_inner_shape);
- }
- static real f_H_gap(real * op, real *ap, real time)
- {
- real not_used_var=0;
- return f_h(op[0], op[1], ap[0], ap[1], not_used_var, &pump_gap_outer_shape);
- }
- /* Functions used to calcuate a new position after a rigid body rotation */
- static void find_app_rigid(real * ap, real * app, real dtheta)
- {
- app[0]=ap[0]*cos(dtheta)-ap[1]*sin(dtheta);
- app[1]=ap[0]*sin(dtheta)+ap[1]*cos(dtheta);
- }
- /* For a give ap, it finds the new position appp.
- According to note page 1, the inputs should be ap, op, dt, functions used to calculate h and H. time is
- also given as an input in case h and H are functions of time. But for this case, h and H are not.
- */
- static void find_appp(real * appp, real * ap, real *op, real time, real dt, real (*ff_h)(real *, real *, real), real (*ff_H)(real *, real *, real))
- {
- real H_t, H_t_dt, h_t, h_t_dt, dtheta, abs_oppappp, abs_opap;
- real app[3], opp[2], opap[2], temp[2];
- app[2] = ap[2];
- dtheta=RPM*M_PI/30*dt;
- find_app_rigid(ap, app, dtheta); /* Find app, which is after a rigid body rotation of dtheta */
- find_app_rigid(op, opp, dtheta); /* Find opp, which is after a rigid body rotation of dtheta */
- h_t = ff_h(op, ap, time);
- h_t_dt= ff_h(opp, app, time+dt);
- H_t = ff_H(op, ap, time);
- H_t_dt= ff_H(opp, app, time+dt);
- abs_opap = sqrt(pow(op[0]-ap[0],2)+pow(op[1]-ap[1],2));
- abs_oppappp=(abs_opap-h_t)/(H_t-h_t)*(H_t_dt-h_t_dt)+h_t_dt;
- opap[0]=ap[0]-op[0];
- opap[1]=ap[1]-op[1];
- temp[0]=opap[0]/abs_opap;
- temp[1]=opap[1]/abs_opap;
- appp[0]=(temp[0]*cos(dtheta)-temp[1]*sin(dtheta))*abs_oppappp+opp[0];
- appp[1]=(temp[0]*sin(dtheta)+temp[1]*cos(dtheta))*abs_oppappp+opp[1];
- }
- DEFINE_GRID_MOTION(vane_pump_gap, domain, dt, time, dtime)
- {
- cell_t c;
- Thread *tc = DT_THREAD ((Dynamic_Thread *)dt);
- int n;
- Node *v;
- real ap[2], op[2], appp[2];
- SET_DEFORMING_THREAD_FLAG (tc);
- begin_c_loop(c, tc)
- {
- c_node_loop(c, tc, n)
- {
- v = C_NODE(c, tc, n);
- if (NODE_POS_NEED_UPDATE(v))
- {
- NODE_POS_UPDATED(v);
- ap[0]=NODE_X(v);
- ap[1]=NODE_Y(v);
- f_Op_gap(op, op+1, ap[0], ap[1], time-dtime, C_UDMI(c, tc, UDM_sector));
- find_appp(appp, ap, op, time-dtime, dtime, f_h_gap, f_H_gap);
- NODE_X(v)=appp[0];
- NODE_Y(v)=appp[1];
- }
- }
- Update_Cell_Metrics (c, tc);
- }
- end_c_loop (c, tc);
- }
- DEFINE_GRID_MOTION(vane_pump_core, domain, dt, time, dtime)
- {
- cell_t c;
- Thread *tc = DT_THREAD ((Dynamic_Thread *)dt);
- int n;
- Node *v;
- real ap[3], op[2], appp[2];
- /* set deforming flags */
- SET_DEFORMING_THREAD_FLAG (tc);
- begin_c_loop(c, tc)
- {
- c_node_loop(c, tc, n)
- {
- v = C_NODE(c, tc, n);
- if (NODE_POS_NEED_UPDATE(v))
- {
- NODE_POS_UPDATED(v);
- ap[0]=NODE_X(v);
- ap[1]=NODE_Y(v);
- ap[2]=NODE_Z(v);
- f_Op_core(op, op+1, time-dtime, C_UDMI(c, tc, UDM_sector));
- find_appp(appp, ap, op, time-dtime, dtime, f_h_core, f_H_core);
- NODE_X(v)=appp[0];
- NODE_Y(v)=appp[1];
- }
- }
- Update_Cell_Metrics (c, tc);
- }
- end_c_loop (c, tc);
- }
- DEFINE_GRID_MOTION(walls, domain, dt, time, dtime)
- {
- }
- static void intialize_one_shape(struct shape * one_shape, enum define_shape shape_type, enum shrink shrink_y_or_n, real * center, real radius, char * file_name)
- {
- one_shape->shape_type=shape_type;
- if (one_shape->shape_type==circle)
- {
- one_shape->center[0]=center[0];
- one_shape->center[1]=center[1];
- one_shape->radius=radius;
- }
- else if (one_shape->shape_type==user_defined)
- {
- FILE *fp_data;
- int i, j;
- one_shape->shrink_y_or_n=shrink_y_or_n;
- #if PARALLEL
- if(I_AM_NODE_HOST_P)
- #endif
- {
- if(!(fp_data=fopen(file_name,"r")))
- {
- Message0("\nCan not open file %s -aborting!!", file_name);
- exit(0);
- }
- fscanf(fp_data, "%d", &(one_shape->number_of_data));
- }
- host_to_node_int_1(one_shape->number_of_data);
- if(!(one_shape->data=(real (*)[2]) calloc(one_shape->number_of_data, 2*sizeof(real))))
- {
- Message0("\nMemory allocatoin failure-aborting!!\n");
- exit(0);
- }
- for(i=0; i<one_shape->number_of_data; i++)
- for(j=0; j<2; j++)
- {
- #if PARALLEL
- if(I_AM_NODE_HOST_P)
- #endif
- {
- #if RP_DOUBLE
- fscanf(fp_data, "%le", &(one_shape->data[i][j]));
- #else
- fscanf(fp_data, "%e", &(one_shape->data[i][j]));
- #endif
- }
- }
- host_to_node_real(&(one_shape->data[0][0]), 2*one_shape->number_of_data);
- #if PARALLEL
- if(I_AM_NODE_HOST_P)
- #endif
- {
- fclose(fp_data);
- }
- }
- else if (one_shape->shape_type==special)
- {
- }
- else
- {
- Message0("\nWrong shape-aborting!!!\n");
- exit(0);
- }
- }
- static void initialize_shape(void)
- {
- real center[2], radius;
- center[0]=-DELTA;
- center[1]=0;
- radius=R;
- intialize_one_shape(&pump_core_outer_shape, OUTER_SHAPE, wo_shrink, center, radius, "data_outer.txt");
- center[0]=0;
- center[1]=0;
- radius=r;
- intialize_one_shape(&pump_core_inner_shape, INNER_SHAPE, wo_shrink, center, radius, "data_inner.txt");
- center[0]=-DELTA;
- center[1]=0;
- radius=R;
- intialize_one_shape(&pump_gap_outer_shape, OUTER_SHAPE, wo_shrink, center, radius, "data_outer.txt");
- center[0]=-DELTA;
- center[1]=0;
- radius=R-GAP;
- intialize_one_shape(&pump_gap_inner_shape, OUTER_SHAPE, w_shrink, center, radius, "data_outer.txt");
- return;
- }
- /* Initialization function is used to save the inital angle alfa used in op calculation.
- It is also needed to read the outer contour for user_defined shape. */
- DEFINE_INIT(init_sector, domain)
- {
- Thread *tc;
- cell_t c;
- Node *v;
- int j;
- real angle;
- /* Initialization can only be done at the beginning */
- Message0("\n\n**************************************************************************");
- Message0("\n\n Warning: you are initializing the flow. For the vane pump application,");
- Message0("\n you can only do this when the mesh is at the initial position. Otherwise,");
- Message0("\n initialization may cause mesh motion failure!!!\n");
- Message0("\n After initialization, please display the UDM-0 using cell value to verify");
- Message0("\n that the sector numbers are calculated correctly!\n");
- Message0("\n**************************************************************************\n");
- /* Initialize the Core sector number*/
- tc=Lookup_Thread(domain, CORE_ID);
- begin_c_loop_int(c, tc)
- {
- v = C_NODE(c, tc, 0);
- j = my_atan2(NODE_Y(v),NODE_X(v))/(2*M_PI/N_VANE);
- C_UDMI(c, tc, UDM_sector) = j;
- }
- end_c_loop_int(c, tc)
- /* Initialize the Gap sector number*/
- if (GAP_ID != 0) {
- tc=Lookup_Thread(domain, GAP_ID);
- begin_c_loop_int(c, tc)
- {
- v = C_NODE(c, tc, 0);
- angle = my_atan2(NODE_Y(v),NODE_X(v));
- if (fabs(angle-2*M_PI)<2*M_PI/N_VANE*0.3)
- angle = 0;
- j = (angle+2*M_PI/N_VANE*0.5)/(2*M_PI/N_VANE);
- C_UDMI(c, tc, UDM_sector) = j;
- }
- end_c_loop_int(c, tc)
- }
- }
- /* print the outer contour for user_defined shape_type*/
- static void print_one_shape(struct shape * one_shape)
- {
- Message0("%d\n", one_shape->shape_type);
- Message0("%d\n", one_shape->shrink_y_or_n);
- if (one_shape->shape_type==circle)
- {
- Message0("%e %e %e\n", one_shape->center[0], one_shape->center[1], one_shape->radius);
- }
- else if(one_shape->shape_type==user_defined)
- {
- Message0("%d\n", one_shape->number_of_data);
- #ifdef DEBUG
- int i, j;
- for(i=0; i<one_shape->number_of_data; i++)
- {
- for(j=0; j<2; j++)
- {
- Message0("%12.4e", one_shape->data[i][j]);
- }
- Message0("\n");
- }
- #endif
- }
- else if(one_shape->shape_type==special)
- {
- Message0("\nSpecial Inner shape\n");
- }
- else
- {
- Message0("\nWrong type-aborting!!\n"); exit(0);
- }
- }
- /* Only node 0 will execute ON_DEMAND udf */
- DEFINE_ON_DEMAND(print_outer_contour)
- {
- Message0("\n******************** Below is the pump core outer infor ********************\n\n");
- print_one_shape(&pump_core_outer_shape);
- Message0("\n******************** Below is the pump core inner infor ********************\n\n");
- print_one_shape(&pump_core_inner_shape);
- if(GAP_ID!=0)
- {
- Message0("\n******************** Below is the pump gap outer infor ********************\n\n");
- print_one_shape(&pump_gap_outer_shape);
- Message0("\n******************** Below is the pump gap inner infor ********************\n\n");
- print_one_shape(&pump_gap_inner_shape);
- }
- }
- /* The input file will be read each time the UDF is loaded */
- DEFINE_EXECUTE_ON_LOADING(on_loading, libudf)
- {
- FILE* fpin;
- Message0("\nON_LOADING udf is executing ...\n");
- #if PARALLEL
- if(I_AM_NODE_HOST_P)
- #endif
- {
- fpin = fopen("input.txt","r");
- if(fpin == NULL)
- Message0("Input file does not exist!");
- }
- #if RP_DOUBLE
- # define REAL_FMT "%lg"
- #else /* RP_DOUBLE */
- # define REAL_FMT "%g"
- #endif /* RP_DOUBLE */
- #if PARALLEL
- if(I_AM_NODE_HOST_P)
- #endif
- {
- fscanf(fpin, "%d",&OUTER_SHAPE);
- fscanf(fpin, "%d",&INNER_SHAPE);
- fscanf(fpin, "%d",&N_VANE);
- fscanf(fpin, "%d %d",&CORE_ID, &GAP_ID);
- fscanf(fpin, REAL_FMT REAL_FMT REAL_FMT, &r, &HALF_VANE_WIDTH, &GAP);
- fscanf(fpin, REAL_FMT REAL_FMT, &R, &DELTA);
- fclose(fpin);
- }
- host_to_node_int_5(OUTER_SHAPE, INNER_SHAPE, N_VANE, CORE_ID, GAP_ID);
- host_to_node_real_5(r, HALF_VANE_WIDTH, GAP, R, DELTA);
- /* Read in any profile data files */
- initialize_shape();
- }
- static void sector_output(int flag, int sector_number)
- {
- char filename[15]="output-";
- char filename_1[1], filename_2[4]=".txt";
- real Operating_Pressure=0;
- real Chamber_Volume_Weighted_Pressure_Sum=0;
- real Chamber_Pressure=0;
- real Chamber_Volume=0;
- FILE *fp;
- Domain *domain;
- Thread *tc;
- cell_t c;
- #if PARALLEL
- if(I_AM_NODE_HOST_P)
- #endif
- {
- filename_1[0]=flag+48;
- strncat(filename, filename_1, 1);
- strcat(filename, ".txt");
- if(N_TIME == 1) {
- fp=fopen(filename,"w+");
- fprintf(fp, " time(s) volume (cc) pressure (Pa) \n");
- fclose(fp);
- }
- fp=fopen(filename, "a");
- }
- domain = Get_Domain(1);
- tc=Lookup_Thread(domain, CORE_ID);
- Operating_Pressure = RP_Get_Real ("operating-pressure");
- begin_c_loop(c,tc)
- {
- if (C_UDMI%
复制代码 |