本案例中要解决的问题如下图所示。
本案例利用UDF定义质量源项与能量源项来模拟膜沸腾现象。计算区域宽0.0389 m,高0.1168 m。介质的饱和温度 K,计算区域底部壁面的温度为510 K。初始时刻,沿Y方向给定一个从Twall到TSAT的线性分布温度,利用下面的函数表达式进行指定:
计算结果如下图所示。
创建两个新材料liquid及vapor,材料属性如下图所示。
image-20220106133257873
UDF中使用了3个UDM,这里需要先设置UDM数量。
案例实际上是投机取巧了的,事实上膜状沸腾模拟最困难的地方是附壁气膜的形成过程,这里直接将气膜以初始化的形式直接给定了。
气相质量源的一般表达形式为:
式中,为通过界面的单位面积的热通量,下标与分别表示气相与液相;为潜热。
这样可以得到气相质量源为:
由于计算域内没有内部质量源,因此液相质量源为:
能量方程的潜热源变为:
界面属性包括表面张力系数0.1 N/m,潜热1e5 J/kg,饱和温度 K。
案例模型的长度尺度是泰勒-罗利不稳定性最危险的波长:
速度尺度为:
因此时间尺度为:
计算域水平宽度为,垂直高度为,网格分辨率为64(水平方向)×192(高度方向)。
汽液界面的初始形状必须受到扰动才能启动气泡的生长。因此,采用一个初始化UDF宏,利用气体填充所有满足以下条件的网格:
式中,为坐标,单位为米。
Nusselt数是表征沸腾换热的一个重要的无量纲值,其定义为:
由于该问题的时间尺度为0.1s,因此时间步长为0.001,即100个时间步长分辨率。总而言之,这个问题应该运行大约1200个时间步才能捕捉到第一个气泡排放。
一些代码注释如下所示。
#include "udf.h"#include "sg.h"#include "sg_mphase.h"#include "flow.h"#include "mem.h"DEFINE_ADJUST(area_density, domain){ Thread *t; Thread **pt; cell_t c; //混合相中循环 mp_thread_loop_c(t, domain, pt) if (FLUID_THREAD_P(t)) { //利用P_PHASE得到主相的Thread指针,P_PHASE宏的值=0 Thread *tp = pt[P_PHASE]; begin_c_loop(c, t) { //计算得到$\namda T \cdot \nabla α$,并将其存储在UDM中方便后面调用 C_UDMI(c, t, 0) = (C_VOF_G(c, tp)[0] * C_T_G(c, t)[0] C_VOF_G(c, tp)[1] * C_T_G(c, t)[1]); } end_c_loop(c, t) }}//气相质量源DEFINE_SOURCE(gas, cell, thread, dS, eqn){ real x[ND_ND]; real source; Thread *tm = THREAD_SUPER_THREAD(thread); Thread **pt = THREAD_SUB_THREADS(tm); //得到k_lα_l k_gα_g real Kl = C_K_L(cell, pt[1]) * C_VOF(cell, pt[1]), Kg = C_K_L(cell, pt[0]) * C_VOF(cell, pt[0]); real L = 1e5; // 根据公式计算得到质量源,这里引用了前面的UDM变量 source = (Kl Kg) * C_UDMI(cell, tm, 0) / L; //将质量源存储在UDM中方便后面调用 C_UDMI(cell, tm, 1) = source; // 顺便计算能量源(等于质量源于潜热的乘积),存储在UDM中后面调用 C_UDMI(cell, tm, 2) = -source * L; // 无法保证负斜率,干脆直接赋零 dS[eqn] = 0; return source;}//液相质量源,利用UDM值进行指定DEFINE_SOURCE(liquid, cell, thread, dS, eqn){ real x[ND_ND]; real source; Thread *tm = THREAD_SUPER_THREAD(thread); Thread **pt = THREAD_SUB_THREADS(tm); source = -C_UDMI(cell, tm, 1); dS[eqn] = 0; return source;}// 计算能量源DEFINE_SOURCE(energy, cell, thread, dS, eqn){ real x[ND_ND]; real source; Thread *tm = thread; source = C_UDMI(cell, tm, 2); dS[eqn] = 0; return source;}//给底部初始化了一层气体区域DEFINE_INIT(my_init_function, domain){ Thread *t; Thread **pt; Thread **st; cell_t c; real xc[ND_ND], y, x; mp_thread_loop_c(t, domain, pt) if (FLUID_THREAD_P(t)) { Thread *tp = pt[P_PHASE]; begin_c_loop(c, t) { C_CENTROID(xc, c, t); x = xc[0]; y = xc[1]; if (y < 0.00292 0.0006 * cos(6.283 * x / 0.0778)) C_VOF(c, tp) = 1; else C_VOF(c, tp) = 0; } end_c_loop(c, t) }}
免责声明:本文系网络转载或改编,未找到原创作者,版权归原作者所有。如涉及版权,请联系删