一、问题背景

最近做的一个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_HOST
    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)
    #endif
}

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

小标题啥意思呢？

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

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

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

        #if !RP_HOST
        sheet_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_HOST
        sheet_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_HOST
    Thread *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_NODE
    Message("UDM1 is initialized successfully.\n");

    #endif
}