给流场以分布的柱坐标形式初场

Rainyboy

较早的一项工作，我已不从事这方面的工作，贴出来供大家参考吧。

在旋转机械的非定常计算中，往往希望指定一个合适的初场，这将有利于计算的收敛、甚至决定最终收敛的结果。

FLUENT的GUI不能提以供柱坐标表示的初始条件，导致初始化的整周叶片通道中，部分通道流场中有不应该出现的涡。

通过UDF可以完美的解决这个问题，代码如下：

1. /************************************************************************/
   
2. /* 初始化流场：给流场以分布的柱坐标形式初场 */
   
3. /* 作者： 范 雨
   
4. /* 创建时间： 2009.9.21
   
5. /* 历次修改：
   
6. /************************************************************************/
   
7. 
   
8. 
   
9. #include "udf.h"
   
10. #include <stdio.h>
    
11. #include <math.h>
    
12. 
    
13. 
    
14. 
    
15. 
    
16. //速度常量
    
17. #define VELOCITY_R_DIRECTION 0.0
    
18. #define VELOCITY_A_DIRECTION -146.0
    
19. #define VELOCITY_Z_DIRECTION 167.0
    
20. 
    
21. //定义圆周率
    
22. #define PI 3.141592653
    
23. 
    
24. //接近零的数
    
25. #define NEAR_ZERO 1e-15
    
26. 
    
27. //记录柱坐标系向量（或者点）的结构体
    
28. struct _Vector_Rot
    
29. {
    
30.     real r;
    
31.     real angle;
    
32.     real z;
    
33. };
    
34. 
    
35. //记录直角坐标系向量（或者点）的结构体
    
36. struct _Vector_XY
    
37. {
    
38.     real x;
    
39.     real y;
    
40.     real z;
    
41. };
    
42. 
    
43. //柱坐标系向量转直角坐标系
    
44. struct _Vector_XY Rot_To_XY(struct _Vector_Rot RotPossition)
    
45. {
    
46.     struct _Vector_XY XYPossition;
    
47.     XYPossition.x = RotPossition.r * cos(RotPossition.angle);
    
48.     XYPossition.y = RotPossition.r * sin(RotPossition.angle);
    
49.     XYPossition.z = RotPossition.z;
    
50.     return XYPossition;
    
51. }
    
52. 
    
53. //直角坐标系向量转柱坐标系
    
54. struct _Vector_Rot XY_To_Rot(struct _Vector_XY XYPossition)
    
55. {
    
56.     struct _Vector_Rot RotPossition;
    
57.     RotPossition.r = pow(pow(XYPossition.x,2)+pow(XYPossition.y,2),0.5);
    
58.     if (RotPossition.r > -1*NEAR_ZERO && RotPossition.r < NEAR_ZERO)
    
59.     {
    
60.         //原点
    
61.         RotPossition.angle = 0;
    
62.     }
    
63.     else
    
64.     {
    
65.         //不是原点
    
66.         if (XYPossition.y < NEAR_ZERO)
    
67.         {
    
68.             //第三、四象限
    
69.             RotPossition.angle = 2* PI - acos(XYPossition.x/RotPossition.r) ;
    
70.         }
    
71.         else
    
72.         {
    
73.             //第一、二象限
    
74.             RotPossition.angle = acos(XYPossition.x/RotPossition.r);
    
75.         }
    
76.     }
    
77.     RotPossition.z = XYPossition.z;
    
78.     return RotPossition;
    
79. }
    
80. 
    
81. //得到与当前向量（直角坐标）在同一Z位置垂直的单位向量（柱坐标）
    
82. struct _Vector_Rot GetXYVector_CZ_CurrentXYVector(struct _Vector_XY CurrentXYVector)
    
83. {
    
84.     struct _Vector_Rot CurrentRotVector = XY_To_Rot(CurrentXYVector);
    
85.     CurrentRotVector.angle += PI/2.0;
    
86.     return CurrentRotVector;
    
87. }
    
88. 
    
89. //得到当前点（直角坐标）的切向速度（直角坐标）
    
90. struct _Vector_XY GetRotVelocity_ByXYPossition(struct _Vector_XY XYPossition)
    
91. {
    
92.     struct _Vector_Rot CurrentRotVector = GetXYVector_CZ_CurrentXYVector(XYPossition);
    
93.     CurrentRotVector.r = VELOCITY_A_DIRECTION;
    
94.     return Rot_To_XY(CurrentRotVector);
    
95. }
    
96. //得到当前点（直角坐标）的径向速度（直角坐标）
    
97. struct _Vector_XY GetRadixVeclocity_ByXYPossition(struct _Vector_XY XYPossition)
    
98. {
    
99.     struct _Vector_Rot CurrentRadixVector = XY_To_Rot(XYPossition);
    
100.     CurrentRadixVector.r = VELOCITY_R_DIRECTION;
     
101.     return Rot_To_XY(CurrentRadixVector);
     
102. }
     
103. 
     
104. //由点在直角坐标系的坐标得到其速度（直角坐标系表示）
     
105. struct _Vector_XY GetXYVelocity_ByXYPossiontion(struct _Vector_XY XYPossition)
     
106. {
     
107.     struct _Vector_XY CurrentRadixVector_inXY = GetRadixVeclocity_ByXYPossition(XYPossition);
     
108.     struct _Vector_XY CurrentRotVector_inXY = GetRotVelocity_ByXYPossition(XYPossition);
     
109.     struct _Vector_XY CurrentTotalVector_inXY;
     
110.     CurrentTotalVector_inXY.x = CurrentRadixVector_inXY.x + CurrentRotVector_inXY.x;
     
111.     CurrentTotalVector_inXY.y = CurrentRadixVector_inXY.y + CurrentRotVector_inXY.y;
     
112.     CurrentTotalVector_inXY.z = VELOCITY_Z_DIRECTION;
     
113.     return CurrentTotalVector_inXY;
     
114. }
     
115. 
     
116. 
     
117. /************************************************************************/
     
118. /* 在Fluent UDF 中的代码 */
     
119. /************************************************************************/
     
120. 
     
121. DEFINE_INIT(Further_Init,domain)
     
122. {
     
123.     FILE *fp;
     
124.     cell_t Nod_cell;
     
125.     Thread *pNod_Thread;
     
126.     real CellPossition[ND_ND];
     
127. 
     
128.     struct _Vector_XY sPossition;
     
129.     struct _Vector_XY sVelocity;
     
130.     fp = fopen("data.txt","w");
     
131.     Message("/************************************************************************/\n");
     
132.     Message("/* 初始化流场：给流场以分布的柱坐标形式初场 */\n");
     
133.     Message("/* 作者： 范 雨 */\n");
     
134.     Message("/* 创建时间： 2009.9.21 */\n");
     
135.     Message("/* 测试编号： 3.4001 */\n");
     
136.     Message("/************************************************************************/\n");
     
137.     thread_loop_c(pNod_Thread,domain)
     
138.     {
     
139.         begin_c_loop_all(Nod_cell,pNod_Thread){
     
140.         //得到坐标（直角）
     
141.         C_CENTROID(CellPossition,Nod_cell,pNod_Thread);
     
142.         fprintf(fp,"Cell Possition %f ,%f ,%f\n",CellPossition[0],CellPossition[1],CellPossition[2]);
     
143.         sPossition.x = CellPossition[0];
     
144.         sPossition.y = CellPossition[1];
     
145.         sPossition.z = CellPossition[2];
     
146.         //得到速度（直角）
     
147.         sVelocity = GetXYVelocity_ByXYPossiontion(sPossition);
     
148.         //设置速度（直角）
     
149.         C_U(Nod_cell,pNod_Thread) = sVelocity.x;
     
150.         C_V(Nod_cell,pNod_Thread) = sVelocity.y;
     
151.         C_W(Nod_cell,pNod_Thread) = sVelocity.z;
     
152.         fprintf(fp,"Velocity %f ,%f ,%f\n\n",sVelocity.x,sVelocity.y,sVelocity.z);
     
153.         }end_c_loop_all(Nod_cell,pNod_Thread)
     
154.     }
     
155.     fclose(fp);
     
156. }

复制代码

详细的流程为（附件即是此图，同Visio打开）：


用一个简单的圆筒流场示例，运行该代码后的速度矢量（为避免显示混乱，只显示了内外侧的速度矢量）为：

Heidiadalheid

高人，下载学习一下

likgo

为何下载不了啊

Heidiadalheid

likgo 发表于 2011-1-18 05:15

                               
登录/注册后可看大图

为何下载不了啊

别用迅雷之类的下载软件

Rainyboy

回复 3 # likgo 的帖子

附件就是那个流程图而已……当时是怕上传的图显示不清楚才作为附件的。