|
阀门流动计算,考虑cavitation效应。
在阀门运动之前,需要判断其上的合力大小,
以判断是否需要开启动网格计算。
用SDOF(或6DOF)进行动网格.
因为需要用到并行计算,
所以参照如下UDF算例编写。
可以求出壁面的合力。
【link1】
(注:这个例子是静压的平均值,
稍加修改即可计算合力)
#include "udf.h" DEFINE_ADJUST(face_av,domain) { /* Variables used by serial, host, node versions */ int surface_thread_id=0; real total_area=0.0; real total_force=0.0; /* "Parallelized" Sections */ #if !RP_HOST /* Compile this section for computing processes only (serial and node) since these variables are not available on the host */ Thread* thread; face_t face; real area[ND_ND]; #endif /* !RP_HOST */ /* Get the value of the thread ID from a user-defined Scheme variable */ #if !RP_NODE /* SERIAL or HOST */ surface_thread_id = RP_Get_Integer("pres_av/thread-id"); Message("\nCalculating on Thread # %d\n",surface_thread_id); #endif /* !RP_NODE */ /* To set up this user Scheme variable in cortex type */ /* (rp-var-define 'pres_av/thread-id 2 'integer #f) */ /* Once set up you can change it to another thread's ID using : */ /* (rpsetvar 'pres_av/thread-id 7) */ /* Send the ID value to all the nodes */ host_to_node_int_1(surface_thread_id); /* Does nothing in serial */ #if RP_NODE Message("\nNode %d is calculating on thread # %d\n",myid, surface_thread_id); #endif /* RP_NODE */ #if !RP_HOST /* SERIAL or NODE */ /* thread is only used on compute processes */ thread = Lookup_Thread(domain,surface_thread_id); //获取壁面face的thread begin_f_loop(face,thread) //对face进行遍历搜索,每段上都求力 /* get the area vector and pressure and increment */ /* the total area and total force values for this node */ if (PRINCIPAL_FACE_P(face,thread)) /* Always TRUE in serial version */ { F_AREA(area,face,thread); //获取face的面矢量area,即面的法线(未经归一化) total_area += NV_MAG(area); total_force += NV_MAG(area)*F_P(face,thread);//面积乘以静压,得到力。这如果乘以面矢量的分量(如area[0]是x分量),就得到分力 } end_f_loop(face,thread)Message("Total Area Before Summing %f\n",total_area); Message("Total Normal Force Before Summing %f\n",total_force); # if RP_NODE /* Perform node synchronized actions here Does nothing in Serial */ total_area = PRF_GRSUM1(total_area); total_force = PRF_GRSUM1(total_force); # endif /* RP_NODE */ #endif /* !RP_HOST */ /* Pass the node's total area and pressure to the Host for averaging */ node_to_host_real_2(total_area,total_force); /* Does nothing in SERIAL */ #if !RP_NODE /* SERIAL or HOST */ Message("Total Area After Summing: %f (m2)\n",total_area); Message("Total Normal Force After Summing %f (N)\n",total_force); Message("Average pressure on Surface %d is %f (Pa)\n", surface_thread_id,(total_force/total_area)); #endif /* !RP_NODE */ } |
但是,
这个例子,
仅适用于单相流动。
对于气液两相流,
计算的不对!
主要是因为,
两相流时,
混合静压=液相静压 X 液相体积分数 + 气相静压 X 气相体积分数
找到各相的体积分数(volume fraction),就可以了。
提取cell和face的体积分数参考下面的例子。
【link2】
# include "udf.h" # define domain_ID 2 DEFINE_ADJUST(adjust_gradient, domain) { Thread *t; cell_t c; face_t f; domain = Get_Domain(domain_ID); /* Fill UDS with the variable. */ thread_loop_c (t,domain) { begin_c_loop (c,t) { C_UDSI(c,t,0) = C_VOF(c,t); //获取单元的volume fraction } end_c_loop (c,t) } thread_loop_f (t,domain) { if (THREAD_STORAGE(t,SV_UDS_I(0))!=NULL) begin_f_loop (f,t) { F_UDSI(f,t,0) = F_VOF(f,t); //获取face的olume fraction } end_f_loop (f,t) } } |
找到了体积分数,
在合力前面乘上就行。
但是,还缺了一步。
怎么区分液相和气相的静压呢?
【link1】实际上没有区分。
下面的例子能够找到液相和气相。
【link3】
/***************************************************************** UDF for initializing phase volume fraction ******************************************************************/ #include "udf.h" /* domain pointer that is passed by INIT function is mixture domain */ DEFINE_INIT(my_init_function, mixture_domain) { int phase_domain_index; cell_t cell; Thread *cell_thread; Domain *subdomain; real xc[ND_ND]; /* loop over all subdomains (phases) in the superdomain (mixture) */ sub_domain_loop(subdomain, mixture_domain, phase_domain_index) //对domain下的各相进行遍历,mixture_domain可以通过Get_Domain(1)获得(默认1对应混合相,2为主相,3为次相) { /* loop if secondary phase */ if (DOMAIN_ID(subdomain) == 3) /* loop over all cell threads in the secondary phase domain */ thread_loop_c (cell_thread,subdomain) { /* loop over all cells in secondary phase cell threads */ begin_c_loop_all (cell,cell_thread) { C_CENTROID(xc,cell,cell_thread); if (sqrt(ND_SUM(pow(xc[0] - 0.5,2.), pow(xc[1] - 0.5,2.), pow(xc[2] - 0.5,2.))) < 0.25) /* set volume fraction to 1 for centroid */ C_VOF(cell,cell_thread) = 1.; else /* otherwise initialize to zero */ C_VOF(cell,cell_thread) = 0.; } end_c_loop_all (cell,cell_thread) } } } |
注意:
直接复制上述UDF编译时,可能出现错误。
原因是上面的空格,不是UTF-8格式,需要全部替换。
我用的小例子如下:
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-4-20 11:14
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社