Fluent UDF 任意物理量梯度计算

一 前言：

使用Fluent UDF进行二次开发时，常常需用到物理量的梯度、二阶导数等，而这些在Fluent中并未直接提供，因此需要用户自己计算。本文在UDF中采用最小二乘法，从基础物理量出发，完成物理量的梯度计算。

下图 给出了UDF计算结果与Fluent直接计算结果的对比结果，算例基于三维溃堤算例。UDF计算结果与Fluent计算结果分布趋势与数值基本一致，细微差异可能在于Fluent使用了一些重构或者限制器。

图 1 Fluent结果与UDF结果对比

二 UDF源代码

以下为源代码（可以在线查看，也可以下载源代码），基于该代码，用户不用再通过UDS的方式，间接获取所需梯度。

/* ******************************************************************

* 文件名称： 3D_Gradient.c    

* 内容摘要：

*         对任意存储于网格中心的物理量计算梯度，兼容二维或三维模型

* 其它说明： 

*       1，采用最小二乘法计算控制体梯度

*       2，采用克莱默法则计算解线性方程组

*       3，测试结果与Fluent的结果存在微小差异，请谨慎使用

 ******************************************************************** */

#include"udf.h"

#include"metric.h"

#define Cell_Len 30

/* 对于网格单元（cell）的每个侧面

*    假定有3个变量记录三个方向距离差值，1个记录距离平方根，

*    1个记录物理量差值，采用（3+1+1）*n的结构记录数据，

*    其中n表征网格单元面数量

*/

DEFINE_ON_DEMAND(My_gradient)

{

  Domain *d;

  Thread *t,*t0,*t1;

  cell_t  c, c0, c1;

  real x0[ND_ND],x1[ND_ND];

  double A[ND_ND][ND_ND],B[ND_ND];

  double temp=0.0,tempx=0.0,tempy=0.0,tempz=0.0;

  real Coef[Cell_Len]={0.0};        

  face_t f;

  int n,i,j;

  d = Get_Domain(1);

  thread_loop_c(t,d)

  {

    begin_c_loop(c,t)

    {

    for(i=0;i<Cell_Len;i++)

    {

Coef[i]=0.0;  

    }

    c_face_loop(c, t, n)                             

    { 

    f = C_FACE(c,t,n);

            t1= C_FACE_THREAD(c,t,n);

    if(THREAD_TYPE(t1)==THREAD_F_INTERIOR)