一、问题背景

最近做的一个fluent仿真算例里用到源项宏，源项宏里用UDM定义了树脂固化度场。

在迭代计算的过程中，UDM的取值发散成了无穷大inf（第一次计算取值是NAN），如下图所示。

由于每一次迭代计算过程中，需要用到以前的数据进行累加，所以某一步计算结果是inf，后面所有结果都会是无穷大。

本文介绍一种方法，在改动导致inf错误的代码的同时，可以不用从头开始计算（前提是已经迭代计算的流动时间占总时间的极小比例）。

二、解决办法

2.1 检查是否有UDF数据结构定义在了HOST节点中

以下是出错的源码，可以看到我原来将面网格数据结构定义在了host节点中，这样很容易出错。

因为host节点不存储任何网格数据，强行计算就可能会导致除以0而发散的悲剧。

notes：以下的代码只是profile宏，对于包括头文件的代码以及定义全局变量bou_temp 的代码被省略（实际存在）。

DEFINE_PROFILE(Inlet_Temp, t, i)
{real time = CURRENT_TIME;face_t f;if(time<=7200){bou_temp = 300 + 0.025*time;}else if(time>7200 && time<=18000){bou_temp = 480;}else if(time>18000 && time<=28800){bou_temp = 480 - (time-18000)/60;}else bou_temp = 300;begin_f_loop(f, t){if (PRINCIPAL_FACE_P(f, t)){F_PROFILE(f, t, i) = bou_temp;}}end_f_loop(f, t)
}

改成如下形式。

DEFINE_PROFILE(Inlet_Temp, t, i)
{#if !RP_HOSTreal time = CURRENT_TIME;face_t f;if(time<=7200){bou_temp = 300 + 0.025*time;}else if(time>7200 && time<=18000){bou_temp = 480;}else if(time>18000 && time<=28800){bou_temp = 480 - (time-18000)/60;}else bou_temp = 300;begin_f_loop(f, t){if (PRINCIPAL_FACE_P(f, t)){F_PROFILE(f, t, i) = bou_temp;}}end_f_loop(f, t)#endif
}

2.2 计算面积或体积权重平均物理量时分母是否偏小

小标题啥意思呢？

例如说一个实体中有3个网格单元，你需要计算这个实体的体积平均温度，公式如下。

但是你在写代码的过程中，可能将分子计算成了V1或V2或V3，而如果这个实体中计算节点足够多，你的实体足够大，单个V占整个总体积的比例越小，你的T_ave就会越趋于发散。

我就犯了如此大错，代码如下：

        #if !RP_HOSTsheet_ct = Lookup_Thread(Get_Domain(1), sheet_id);begin_c_loop_int(c, sheet_ct){curing_ratio_weighted_sum += C_VOLUME(c, sheet_ct) * C_UDMI(c, sheet_ct, 0);heat_weighted_sum += C_VOLUME(c, sheet_ct) * C_UDMI(c, sheet_ct, 3);volume_sum += C_VOLUME(c, sheet_ct);}end_c_loop_int(c, sheet_ct)curing_ratio_weighted_sum = PRF_GRSUM1(curing_ratio_weighted_sum);heat_weighted_sum = PRF_GRSUM1(heat_weighted_sum);curing_ratio_ave = curing_ratio_weighted_sum / volume_sum;heat_ave = heat_weighted_sum / volume_sum;End_iter_time = CURRENT_TIME;#endif

我没有将volume_sum用全局约简宏将所有计算节点上的总体积求和起来，因而就容易趋于发散。

代码改成下面的形式。

        #if !RP_HOSTsheet_ct = Lookup_Thread(Get_Domain(1), sheet_id);begin_c_loop_int(c, sheet_ct){curing_ratio_weighted_sum += C_VOLUME(c, sheet_ct) * C_UDMI(c, sheet_ct, 0);heat_weighted_sum += C_VOLUME(c, sheet_ct) * C_UDMI(c, sheet_ct, 3);volume_sum += C_VOLUME(c, sheet_ct);}end_c_loop_int(c, sheet_ct)curing_ratio_weighted_sum = PRF_GRSUM1(curing_ratio_weighted_sum);heat_weighted_sum = PRF_GRSUM1(heat_weighted_sum);volume_sum = PRF_GRSUM1(volume_sum);curing_ratio_ave = curing_ratio_weighted_sum / volume_sum;heat_ave = heat_weighted_sum / volume_sum;End_iter_time = CURRENT_TIME;#endif

三、如何无需从头开始计算就转inf为正常值

关于curing_ratio这个UDM变量，我在源项宏中对它的迭代计算公式是curing_ratio=curing_ratio+curing_ratio_rate*timestep，很显然如果curing_ratio有一步算出来是inf发散的，那么后面的所有迭代步就都是发散的。

因为已经变成inf（并且第二步就已经发散了），我不知道它迭代到当前步是什么值，所以最终就直接决定将其赋值为0，毕竟当前步和初始步之间的步数占总步数的千分之一还不到。

用什么赋值为0呢？

当然是帅气且灵活的DEFINE_ON_DEMAND宏啦，代码如下。

DEFINE_ON_DEMAND(init_UDM1)
{#if !RP_HOSTThread *sheet_ct;int sheet_id = sheet_zone_id;cell_t c;sheet_ct = Lookup_Thread(Get_Domain(1), sheet_id);begin_c_loop_int(c, sheet_ct){C_UDMI(c, sheet_ct, 1) = 0.0;}end_c_loop_int(c, sheet_ct)#endif#if !RP_NODEMessage("UDM1 is initialized successfully.\\n");#endif
}