给流场以分布的柱坐标形式初场
本帖最后由 Rainyboy 于 2010-10-29 11:07 编辑较早的一项工作，我已不从事这方面的工作，贴出来供大家参考吧。
在旋转机械的非定常计算中，往往希望指定一个合适的初场，这将有利于计算的收敛、甚至决定最终收敛的结果。
FLUENT的GUI不能提以供柱坐标表示的初始条件，导致初始化的整周叶片通道中，部分通道流场中有不应该出现的涡。
通过UDF可以完美的解决这个问题，代码如下：
/************************************************************************/
/* 初始化流场：给流场以分布的柱坐标形式初场 */
/* 作者： 范 雨
/* 创建时间： 2009.9.21
/* 历次修改：
/************************************************************************/
#include "udf.h"
#include <stdio.h>
#include <math.h>
//速度常量
#define VELOCITY_R_DIRECTION 0.0
#define VELOCITY_A_DIRECTION -146.0
#define VELOCITY_Z_DIRECTION 167.0
//定义圆周率
#define PI 3.141592653
//接近零的数
#define NEAR_ZERO 1e-15
//记录柱坐标系向量（或者点）的结构体
struct _Vector_Rot
{
real r;
real angle;
real z;
};
//记录直角坐标系向量（或者点）的结构体
struct _Vector_XY
{
real x;
real y;
real z;
};
//柱坐标系向量转直角坐标系
struct _Vector_XY Rot_To_XY(struct _Vector_Rot RotPossition)
{
struct _Vector_XY XYPossition;
XYPossition.x = RotPossition.r * cos(RotPossition.angle);
XYPossition.y = RotPossition.r * sin(RotPossition.angle);
XYPossition.z = RotPossition.z;
return XYPossition;
}
//直角坐标系向量转柱坐标系
struct _Vector_Rot XY_To_Rot(struct _Vector_XY XYPossition)
{
struct _Vector_Rot RotPossition;
RotPossition.r = pow(pow(XYPossition.x,2)+pow(XYPossition.y,2),0.5);
if (RotPossition.r > -1*NEAR_ZERO && RotPossition.r < NEAR_ZERO)
{
//原点
RotPossition.angle = 0;
}
else
{
//不是原点
if (XYPossition.y < NEAR_ZERO)
{
//第三、四象限
RotPossition.angle = 2* PI - acos(XYPossition.x/RotPossition.r) ;
}
else
{
//第一、二象限
RotPossition.angle = acos(XYPossition.x/RotPossition.r);
}
}
RotPossition.z = XYPossition.z;
return RotPossition;
}
//得到与当前向量（直角坐标）在同一Z位置垂直的单位向量（柱坐标）
struct _Vector_Rot GetXYVector_CZ_CurrentXYVector(struct _Vector_XY CurrentXYVector)
{
struct _Vector_Rot CurrentRotVector = XY_To_Rot(CurrentXYVector);
CurrentRotVector.angle += PI/2.0;
return CurrentRotVector;
}
//得到当前点（直角坐标）的切向速度（直角坐标）
struct _Vector_XY GetRotVelocity_ByXYPossition(struct _Vector_XY XYPossition)
{
struct _Vector_Rot CurrentRotVector = GetXYVector_CZ_CurrentXYVector(XYPossition);
CurrentRotVector.r = VELOCITY_A_DIRECTION;
return Rot_To_XY(CurrentRotVector);
}
//得到当前点（直角坐标）的径向速度（直角坐标）
struct _Vector_XY GetRadixVeclocity_ByXYPossition(struct _Vector_XY XYPossition)
{
struct _Vector_Rot CurrentRadixVector = XY_To_Rot(XYPossition);
CurrentRadixVector.r = VELOCITY_R_DIRECTION;
return Rot_To_XY(CurrentRadixVector);
}
//由点在直角坐标系的坐标得到其速度（直角坐标系表示）
struct _Vector_XY GetXYVelocity_ByXYPossiontion(struct _Vector_XY XYPossition)
{
struct _Vector_XY CurrentRadixVector_inXY = GetRadixVeclocity_ByXYPossition(XYPossition);
struct _Vector_XY CurrentRotVector_inXY = GetRotVelocity_ByXYPossition(XYPossition);
struct _Vector_XY CurrentTotalVector_inXY;
CurrentTotalVector_inXY.x = CurrentRadixVector_inXY.x + CurrentRotVector_inXY.x;
CurrentTotalVector_inXY.y = CurrentRadixVector_inXY.y + CurrentRotVector_inXY.y;
CurrentTotalVector_inXY.z = VELOCITY_Z_DIRECTION;
return CurrentTotalVector_inXY;
}
/************************************************************************/
/* 在Fluent UDF 中的代码 */
/************************************************************************/
DEFINE_INIT(Further_Init,domain)
{
FILE *fp;
cell_t Nod_cell;
Thread *pNod_Thread;
real CellPossition;
struct _Vector_XY sPossition;
struct _Vector_XY sVelocity;
fp = fopen("data.txt","w");
Message("/************************************************************************/\n");
Message("/* 初始化流场：给流场以分布的柱坐标形式初场 */\n");
Message("/* 作者： 范 雨 */\n");
Message("/* 创建时间： 2009.9.21 */\n");
Message("/* 测试编号： 3.4001 */\n");
Message("/************************************************************************/\n");
thread_loop_c(pNod_Thread,domain)
{
begin_c_loop_all(Nod_cell,pNod_Thread){
//得到坐标（直角）
C_CENTROID(CellPossition,Nod_cell,pNod_Thread);
fprintf(fp,"Cell Possition %f ,%f ,%f\n",CellPossition,CellPossition,CellPossition);
sPossition.x = CellPossition;
sPossition.y = CellPossition;
sPossition.z = CellPossition;
//得到速度（直角）
sVelocity = GetXYVelocity_ByXYPossiontion(sPossition);
//设置速度（直角）
C_U(Nod_cell,pNod_Thread) = sVelocity.x;
C_V(Nod_cell,pNod_Thread) = sVelocity.y;
C_W(Nod_cell,pNod_Thread) = sVelocity.z;
fprintf(fp,"Velocity %f ,%f ,%f\n\n",sVelocity.x,sVelocity.y,sVelocity.z);
}end_c_loop_all(Nod_cell,pNod_Thread)
}
fclose(fp);
}
详细的流程为（附件即是此图，同Visio打开）：
用一个简单的圆筒流场示例，运行该代码后的速度矢量（为避免显示混乱，只显示了内外侧的速度矢量）为：
高人，下载学习一下 为何下载不了啊 likgo 发表于 2011-1-18 05:15 static/image/common/back.gif
为何下载不了啊
别用迅雷之类的下载软件 回复 3 # likgo 的帖子
附件就是那个流程图而已……当时是怕上传的图显示不清楚才作为附件的。
页:
[1]