基于Fluent中DPM的水滴蒸发冷凝和沸腾过程中传热传质的规律和程序实现

基于Fluent中DPM的水滴蒸发、冷凝和沸腾过程中传热传
质的规律和程序实现
0 引言
本文将详细讲述基于Fluent软件中DPM的水滴蒸发、冷凝和沸腾的控制方程和程序实现。

在包含水滴的多相流中（离散相和连续相）水滴（离散相）会和湿空气（连续相）发生传热传质。

对于传热过程，包括湿空气和水滴表面的对流换热和水滴蒸发和沸腾时的相变热。

对于传质过程，包括湿空气中水蒸气在液滴上冷凝，液滴的蒸发和沸腾。

请注意，在本文中，湿空气中的水只能在水滴上冷凝，不能在没有水滴的情况下生成水滴，且传热过程中不考虑辐射换热并假设液滴的温度试均匀的。

本文最后提供的程序仅适用于ANSYS19.1及其兼容版本。

为了统一词汇，湿空气（连续相）中的组成成分将表示为“组分(species)”，液滴（离散相）中的组成成分将表示为“成分(component)”。

1 传热
液滴和湿空气的传热主要体现在两方面，一方面是液滴与湿空气之间的对流换热，另一方面是液滴和湿空气之间发生传质时的换热。

下式中，等号右侧第一项为对流换热热量，第二项为传质换热量。

m p c p dT p
dt
=ℎA p(T∞−T p)−
dm p
dt
ℎfg
1.1对流换热
湿空气和液滴之间的对流换热遵循牛顿冷却定律。

液滴与湿空气之间的对流换热流量用下式表示：
m p c p dT p
dt
=ℎA p(T∞−T p)
式中：
m p：液滴（离散相）的质量，kg;
c p：液滴的比热容，J/kg∙K−1;
T p：滴液的温度，K；
ℎ：对流换热系数，W/m2∙K−1;
A p：颗粒表面积，A p=πd p2，m2；
d p：颗粒直径，m；
T∞：湿空气的温度，K。

上式中对流换热系数由湿空气和液滴之间的努塞尔数计算得到，液滴和湿空气之间的努塞尔数由下式给出：
Nu=ℎd p
k
=2+0.6×Re
1
2×Pr
1
3
式中：
k：湿空气的导热系数，W/m∙K−1；
Re：雷诺数；
Pr：湿空气的普朗特数，Pr=μc p
k。

其中雷诺数由下式计算：
Re=ρd p|u p⃗⃗⃗⃗ −u∞
⃗⃗⃗⃗⃗ |
μ
式中：
ρ：湿空气的密度，kg/m3；
u p⃗⃗⃗⃗ ：液滴的速度，m/s；
u∞
⃗⃗⃗⃗⃗ ：湿空气的速度，m/s；
μ：湿空气的动力粘度kg/m∙s。

1.2传热换热
传质换热由液滴的质量变化率和相变焓（汽化潜热）计算。

m p c p dT p
dt
=−
dm p
dt
ℎfg
式中：
ℎfg：相变焓，J/kg。

2 传质
液滴与湿空气之间的传质由三部分构成，蒸发、冷凝和沸腾，其中蒸发和冷凝由扩散定律控制，沸腾时则所有吸收的热量用于相变。

2.1 蒸发与冷凝
当液滴的温度位于沸点以下或者液滴表面水的蒸汽压小于环境压力时，液滴不沸腾，只发生蒸发和冷凝。

T vap≤T p<T bp
液滴发生蒸发或冷凝的方向，或者说水相变的方向由液滴表面和湿空气中水分子的物质的量浓度控制，即湿空气中水分子的物质的量浓度小于液滴表面水分子物质的量浓度时发生蒸发，反之发生冷凝。

当蒸发和冷凝的传质速率不大时，传质速率由扩散定律控制。

液滴的蒸发和湿空气中水的冷凝的传质速率由下式给出：
dm p
dt
=N i A p M w,i
式中：
N i：传质通量，kmol/m2∙s−1；
M w,i：第i组分的摩尔质量，kg/kmol。

上式中传质通量N i决定了传质速率，传质通量由下式给出：
N i=k c(C i,s−C i,∞)
式中：
k c：传质系数，m/s；
C i,s：液滴表面的水蒸汽浓度，kmol/m3；
C i,∞：湿空气的水蒸汽浓度，kmol/m3。

其中，液滴表面的水蒸气浓度即为液滴温度下的水蒸汽饱和浓度。

液滴表面的饱和水蒸汽浓度由理想气体状态方程给出。

C i,s=p sat(T p) RT p
式中：
p sat(T p)：处于液滴温度T p下的水饱和蒸汽压，kPa；
R：通用气体常数，8.314m3Pa/mol∙K−1。

液滴所在地的湿空气中水的物质的量浓度同样由理想气体状态方程给出。

C i,∞=X i
p RT∞
式中：
X i：液滴所在地湿空气中水的摩尔分数；
p：液滴所在地湿空气的绝对压力，kPa。

传质系数k c则与舍伍德数有下式的关系。

SℎAB=k c d p
D i,m
=2+0.6Re
1
2Sc
1
3
式中：
SℎAB：舍伍德数；
D i,m：湿空气中水的扩散系数，m2/s；
Sc：施密特数。

其中施密特数由下式计算：
Sc=
μρD i,m
式中：
μ：湿空气的动力粘度，kg/m∙s；
ρ：湿空气的密度，kg/m3。

当传质速率较大时，传质由对流/扩散控制。

由扩散/对流控制的液滴传质速率由下式给出：
dm p
dt
=k c A pρln(1+B m)
式中B m为散裂质量数，散裂质量数用下式计算。

B m=Y i,s−Y i,∞1−Y i,s
式中：
Y i,s：液滴表面的水蒸汽质量分数；
Y i,∞：湿空气中水的质量分数。

2.2沸腾
当液滴的温度达到沸点T bp，或者滴液表面水的蒸汽压大于环境压力时，液滴发生沸腾，此时液滴由湿空气中得到的热量将全部用于相变。

T p≥T bp
沸腾传质速率由下式给出：
dm p
dt
ℎfg= ℎA p(T∞−T p)
3 程序实现
下列程序是用ANSYS19.1的UDF Manual文档中DEFINE_DPM_HEAT_MASS宏的示例程序改写的，源程序请参见相关文档。

该宏仅可以挂载在DPM模型的multicomponent颗粒中，具体的挂载方法请参见相关文档。

请注意，挂载UDF时一定要保持injection中的LAWS项打开。

当UDF 挂载后，原来的multicomponent-LAW虽然不起作用，但是要保持打开。

在创建湿空气时，要把water-vapor组分放在从上往下数第一个。

/***********************************************************************
UDF for defining the heat and mass transport for
multicomponent particle vaporization
水要放在气相第一位
***********************************************************************/
#include"udf.h"
#include"dpm_mem.h"
DEFINE_DPM_HEAT_MASS(multivap, tp, Cp, hgas, hvap, cvap_surf, Z, dydt, dzdt)
{
int ns;
Material *sp;
real P_total = 0.0; /* vapor pressure */
Thread *t0 = TP_CELL_THREAD(tp); /* thread where the particle is in*/
Material *gas_mix = THREAD_MATERIAL(DPM_THREAD(t0, tp)); /* gas mixture material */ Material *cond_mix = TP_MATERIAL(tp); /* particle mixture material*/
cphase_state_t *c = &(tp->cphase[0]); /* cell information of particle location*/ real molwt[MAX_SPE_EQNS]; /* molecular weight of gas species */
real Tp = TP_T(tp); /* particle temperature K*/
real mp = TP_MASS(tp); /* particle mass g*/
real molwt_bulk = 0.; /* average molecular weight in bulk gas */
real Dp = DPM_DIAM_FROM_VOL(mp / TP_RHO(tp)); /* particle diameter m*/
real Ap = DPM_AREA(Dp); /* particle surface m^3*/
real Pr = c->sHeat * c->mu / c->tCond; /* Prandtl number */
real Nu = 2.0 + 0.6 * sqrt(tp->Re) * pow(Pr, 1. / 3.); /* Nusselt number */
real h = Nu * c->tCond / Dp; /* Heat transfer coefficient*/
real dh_dt = h * (c->temp - Tp) * Ap; /* heat source term*/
dydt[0] += dh_dt / (mp * Cp);
dzdt->energy -= dh_dt;
mixture_species_loop(gas_mix, sp, ns)
{
molwt[ns] = MATERIAL_PROP(sp, PROP_mwi); /* molecular weight of gas species */
molwt_bulk += c->yi[ns] / molwt[ns]; /* 平均摩尔质量的倒数*/
}
/* prevent division by zero */
molwt_bulk = MAX(molwt_bulk, DPM_SMALL);
/*蒸发和冷凝*/
int gas_index = TP_COMPONENT_INDEX_I(tp, 0); /* gas species index of vaporization */
/* condensed material */
Material * cond_c = MIXTURE_COMPONENT(cond_mix, 0);
/* diffusion coefficient m^2/s */
real D = DPM_BINARY_DIFFUSIVITY(tp, cond_c, TP_T(tp));
/* Schmidt number 施密特数，液滴的扩散系数D 湿空气的粘度kg/ms 湿空气的密度kg/m^3*/ real Sc = c->mu / (c->rho * D);
/* mass transfer coefficient m/s*/
real k = (2. + 0.6 * sqrt(tp->Re) * pow(Sc, 1. / 3.)) * D / Dp;
/* bulk gas concentration (ideal gas) kmol/m^3*/
real cvap_bulk = c->pressure / UNIVERSAL_GAS_CONSTANT / c->temp * c->yi[gas_index] / molwt_bulk / solver_par.molWeight[gas_index];
/* vaporization rate 传质速率kg/s*/
real vap_rate = k * molwt[gas_index] * Ap * (cvap_surf[0] - cvap_bulk);
/*颗粒质量变化率kg/s*/
dydt[1] -= vap_rate;
/*气相质量变化率kg/s*/
dzdt->species[gas_index] += vap_rate;
/* dT/dt = dh/dt / (m Cp)颗粒温度变化率K/s*/
dydt[0] -= hvap[gas_index] * vap_rate / (mp * Cp);
/* gas enthalpy source term 气相能量变化率J/s*/
dzdt->energy += hgas[gas_index] * vap_rate;
/*water boiling*/
P_total = cvap_surf[0] * Z * UNIVERSAL_GAS_CONSTANT * Tp;/*颗粒表面蒸汽压*/
if (P_total > c->pressure && dydt[0] > 0.)/*当表面蒸汽压大于颗粒所在位置总压且仍吸热则沸腾*/
{
/*沸腾吸热热量J*/
real h_boil = dydt[0] * mp * Cp;
/* keep particle temperature constant */
dydt[0] = 0.;
int gas_index = TP_COMPONENT_INDEX_I(tp, 0);
/*沸腾蒸发率，吸热量全部转化为为相变热*/
real boil_rate = h_boil / hvap[gas_index];
/* particle component mass source term */
dydt[1] -= boil_rate;
/* fluid species source */
dzdt->species[gas_index] += boil_rate;
/* fluid energy source */
dzdt->energy += hgas[gas_index] * boil_rate;
}
}。