求助 泵动网格的建立
我在FLUENT v6.1动网格教程中发现的一个例子：详见附件看样子非常像我正在做毕设的旋片泵的模型：
http://www.vacuum-guide.com/images/hw_rotary-vane_pneumofore.gif
非常恳求各位大哥大姐，能告诉小弟这个例子具体的动网格设置方法 这个问题温正等人的书上有。 把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.02750 - 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 (xh@fluent.com), and Laz Foley (lnf@fluent.com)
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;
real radius;
real (*data);
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), int number_of_data, real xc, real yc, real xd, real yd,
real *x, real *y)
{
int i, flag;
real pt1, pt2, t1, t2;
for(i=0; i<number_of_data; i++)
{
if(i<number_of_data-1)
{
pt1 = data;
pt1 = data;
pt2 = data;
pt2 = data;
}
else
{
pt1 = data;
pt1 = data;
pt2 = data;
pt2 = data;
}
find_t(pt1, pt1, pt2, pt2, 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, pt1, pt2, pt2, xc, yc, xd, yd, t1, t2, flag); */
if(((flag==1)||(flag==2))&&(t1<=1)&&(t1>=0)&&(t2>=0))
{
(*x) = (1-t1)*pt1 + t1*pt2;
(*y) = (1-t1)*pt1 + t1*pt2;
/* Message0("\n%7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f \n", pt1, pt1, pt2, pt2, 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;
y0 = contour_shape->center;
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, op, ap, ap, ap, &pump_core_inner_shape);
}
static real f_H_core(real * op, real *ap, real time)
{
real not_used_var=0;
return f_h(op, op, ap, ap, 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, op, ap, ap, 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, op, ap, ap, 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=ap*cos(dtheta)-ap*sin(dtheta);
app=ap*sin(dtheta)+ap*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, opp, opap, temp;
app = ap;
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-ap,2)+pow(op-ap,2));
abs_oppappp=(abs_opap-h_t)/(H_t-h_t)*(H_t_dt-h_t_dt)+h_t_dt;
opap=ap-op;
opap=ap-op;
temp=opap/abs_opap;
temp=opap/abs_opap;
appp=(temp*cos(dtheta)-temp*sin(dtheta))*abs_oppappp+opp;
appp=(temp*sin(dtheta)+temp*cos(dtheta))*abs_oppappp+opp;
}
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, op, appp;
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=NODE_X(v);
ap=NODE_Y(v);
f_Op_gap(op, op+1, ap, ap, 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;
NODE_Y(v)=appp;
}
}
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, op, appp;
/* 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=NODE_X(v);
ap=NODE_Y(v);
ap=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;
NODE_Y(v)=appp;
}
}
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=center;
one_shape->center=center;
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 (*)) 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));
#else
fscanf(fp_data, "%e", &(one_shape->data));
#endif
}
}
host_to_node_real(&(one_shape->data), 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, radius;
center=-DELTA;
center=0;
radius=R;
intialize_one_shape(&pump_core_outer_shape, OUTER_SHAPE, wo_shrink, center, radius, "data_outer.txt");
center=0;
center=0;
radius=r;
intialize_one_shape(&pump_core_inner_shape, INNER_SHAPE, wo_shrink, center, radius, "data_inner.txt");
center=-DELTA;
center=0;
radius=R;
intialize_one_shape(&pump_gap_outer_shape, OUTER_SHAPE, wo_shrink, center, radius, "data_outer.txt");
center=-DELTA;
center=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, one_shape->center, 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);
}
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="output-";
char filename_1, filename_2=".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=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% 该同学可能已经毕业了吧。
看到回帖的udf，看得真累。
其实这个动网格不难，思路理清楚之后就慢慢做就是了。
先建立物理模型，搞清楚哪些边界是运动的，是滑动还是转动。
将这些运动的部分建成六面体网格，节点编号要有规律，这样进行动网格策略会很容易操作。
另外，对于此模拟。还有一个比较关键的是高压向低压区泄漏问题，真实情况下是否有其他的密封，比如油封，泄漏速度太大也会导致计算不稳定。
其实动网格这种对软件要求较高的计算，最好用starcd，比较好做。 看到回帖的udf，看得真累。
UDF语言这么复杂？ 有一本书有关于螺旋桨和离心泵的流场仿真：将叶轮及其周围的区域分为一个体（可设定为动），其他区域为另外一个体，两个体之间的面设为interface
这样应该不用编程就可以了。你试试 似乎是 FLUENT高级应用
页:
[1]