DPM添加剪切应力

83136715

想用UDF添加一个由于wall shear stress而产生的力，写了下面的代码，出现错误：line28: subscripted expression is not an array or pointer : double
，请大家帮忙看看。
#include "udf.h"
#include "dpm.h"
#include "storage.h"
#define KVISO 1.48e-06   /*kinematic viscosity*/
#define DP 10.e-06       /*particle diameter*/
#define B 2.e-19         /*B*/
#define Z .4e-9          /*Z*/
#define C .833e-01       /*constant*/

DEFINE_DPM_BODY_FORCE(lift_van, p, i)
{
real lift_force, wall_shear_force, van_force;
real TAU, USTAR;
real lift, adhension;
real area;
real A;
int zone_ID=8;
Thread *t;
face_t f;
cell_t c;
Domain *d;

d=Get_Domain(1);

t=Lookup_Thread(d, zone_ID);
begin_f_loop(f, t)
{
  F_AREA(A,f,t);                             /*出现错误*/
  area=NV_MAG(A);
  wall_shear_force=NV_MAG(F_STORAGE_R_N3V(f, t, SV_WALL_SHEAR));
  TAU=wall_shear_force/area;
  USTAR=sqrt(TAU/C_R(c,t));
}
end_f_loop(f,t)

lift=C_R(c,t)*pow(KVISO,2)*20.90*pow(0.5*USTAR/KVISO,2.31);
adhension=B*DP*C*pow(Z,2);

if(lift>adhension)
  {lift_force=lift;
   lift=0;
   van_force=adhension;
   adhension=0;
  }
else
  continue;        
return((lift_force-van_force)/P_MASS(p));
}

83136715

没人帮忙，呜呜

fty0083

程序本体没法具体去推敲了，大面上看你的A变量定义有问题，fluent里面的面积A是个矢量2，应该定义A[ND_ND]。不过，使用F_AREA(A,f,t)和NV_MAG(A)是没错的。

此外，你程序本身在使用cell的时候有问题，DPM里面可以找到颗粒所在位置的cell和thread。用P_CELL_THREAD(p)表示线，P_CELL(p)表示cell。如此去给定c和t。

那么从你程序的表达来看，是不应该用那个begin_f_loop的。因为，你这里也没有给定f，只是定义了，而没有取值。如果，你想说一个cell的边界面上的话，你也应该先找到cell，再在cell上做界面face循环。也就是说先做上面的给定c和t。具体的计算方法，因为不清楚你的计算本意，因此无法再深究了。

最后，说一下对你这个问题的意见。你想加一个壁面剪切力所产生的力。那么按照你的意思，这个力应该在边界层内才存在。但是，体积力的定义对于所有的颗粒都要作用的，也就是说即使在发展段的颗粒也要计算这个附加力，我不清楚你这里会不会有什么影响。或者，会不会出现在发展段你调用F_STORAGE_R_N3V(f, t, SV_WALL_SHEAR))为空。所以，你再考虑一下具体的code表达。可以加一个判别条件，判别颗粒的位置，再进行计算。比如，对于边界来说，是有zone_id的，fluent帮助也写了一些判别壁面的方法，好好看一看udf manual。

体积力的计算就是返回一个受力大小就可以了，那么这个过程你要先自己给出受力公式，再去求解。那么DPM的这个力来说和定义连续相边界条件是不同的，除非你必须要用到一些界面上的值，否则你是不需要做那种begin_f_loop循环的。先看看自己的模型吧。code重新写一下，那个Domain应该是不需要的。你这里唯一需要注意的应该就是那个壁面剪切应力的取值，具体的取法再思考一下，其他的值都可以直接找到。

83136715

非常感谢你的回复。补充一下我这个code的大体思路：位于viscous sublayer的微粒由于受到wall shear stress的影响而产生了一个lift-off力,如果这个力大于范德华力（van），那么粒子就会发生脱附而进入到发展段，在粒子发生脱附的下一个时刻，这个lift-off力就会消失（范德华力也会消失），不再计算这个受力，也就是像你说的，只考虑边界层。所以在程序段要判断lift-off和van-force的大小比较。
我先按照你给你的意见，改一下我的程序，我也觉得我对udf的学习不够。然后再向你请教。

83136715

对于壁面的判断，我没找到一个简单合适的方法，我本想用lookup_thread(d, zoneID), 但是你说domain的定义应该是不需要的。所以...对于壁面的判断，能不能给点具体的建议。。。thx.

fty0083

关于壁面的判别其实没什么特别好的。从概念上看，对于一个cell，如果包含一个face，这个face在边界上，那么我们就可以说这个cell是边界cell。至于这种边界cell和内部cell的区别这里不去叙述。那么对于你的问题，如果一个cell在边界层内部，那你就应该考虑那个剪切力。一种方法是通过zoneID的方法来判别边界，但是，这个方法实际是欠妥的，因为，边界层内的网格不一定只有边界网格，可能还有好几层网格也在边界层内，所以用zoneID的方法只能确定是边界网格，但是你说明不了边界层。同时，另一种方法也很直接，就是利用点到壁面的距离进行判别，你可以假定一个边界层厚度，如x，那你可以计算当前颗粒位置到壁面的距离，如果小于x，按就应该是在边界层里面了，不过，这种方法对于边界层厚度的取法需要进一步考虑。最后，还有一种方法，你不必拘泥于边界层，所谓壁面剪切力，实际就是因为边界层内流体粘性较大造成的，那么对于内部流体，基本上是没有壁面剪切力的，如此，你可以给定一个连续的剪切力计算公式，作为附加力，这个力越接近壁面，越大，如此，你就不必考虑边界层厚度了。这只是我的几个建议，你还要自己再进一步考虑。

对于计算参数的获取，通过之前我推荐的两个宏，你是可以得到当前particle所在的cell和thread的，那么你就已经拥有了所有当前单元格的流体数据了，具体表达你可以看看cell的相关宏。

此外，如果你想用lookup_thread这个宏的话，给定domain是有意义的。但是先前，你给的程序里面，那个循环的使用本身是有问题的，即face没有给定，所以我觉得你用domain是不必要的。

最后，我建议你不要拘泥于判别壁面，找边界层这些问题，这个剪切力本身就是一直存在的，只不过在边界层里面我们不能忽略而已，所以，最好你能给出一个一般表达，这样对问题也更有说服力。

[ 本帖最后由 fty0083 于 2013-2-25 15:37 编辑 ]

sunning20112

看来大师对DPM中的UDF很了解啊，想请教一下颗粒碰撞壁面时的磨损量怎么用UDF编程呢？要用到碰撞时的颗粒速度、碰撞角、直径，多谢了！

gearboy

整个程序逻辑结构都不正确。随便改了一下，编译通过，但没计算过。整个程序流程应该是：先判断颗粒所在的网格，然后判断该网格是不是壁面的邻接网格，如果是就加你要的力，不是的话就设置为零。

#include "udf.h"

#define KVISO 1.48e-06   /*kinematic viscosity*/
#define DP 10.e-06       /*particle diameter*/
#define BB 2.e-19         /*B*/
#define ZZ .4e-9          /*Z*/
#define CC .833e-01       /*constant*/
const int zone_ID=8;

DEFINE_DPM_BODY_FORCE(lift_van,p,i)
{
        int n;
        real dens,lift_force, wall_shear_force, van_force,TAU, USTAR,lift, adhension,area,result;
        real A[ND_ND];
        face_t wall_f;
        cell_t c;
        Thread *tc,*tf,*wall_tf;       
        Domain*d;
        cxboolean IsCellNeighborToWall=FALSE;

        d=Get_Domain(1);

        //get which cell the particle located in
        c=RP_CELL(&p->cCell);
        tc= RP_THREAD(&p->cCell);
        //judge whether cell is neighbor to wall
        c_face_loop(c,tc,n)
        {
                tf=C_FACE_THREAD(c,tc,n);
                if(zone_ID==THREAD_ID(tf))
                {
                        wall_f=C_FACE(c,tc,n);
                        wall_tf=tf;
                        IsCellNeighborToWall=TRUE;
                }
        }

        if(IsCellNeighborToWall) //if cell is neighbor to wall, then add force
        {
                F_AREA(A,wall_f,wall_tf);
                area=NV_MAG(A);
                wall_shear_force=NV_MAG(F_STORAGE_R_N3V(wall_f,wall_tf,SV_WALL_SHEAR));
                TAU=wall_shear_force/area;
                dens=F_R(wall_f,wall_tf);
                USTAR=sqrt(TAU/dens);
                lift=dens*pow(KVISO,2)*20.90*pow(0.5*USTAR/KVISO,2.31);
                adhension=BB*DP*CC*pow(ZZ,2);
                if(lift>adhension)
                {
                        lift_force=lift;
                        van_force=adhension;
                        result=(lift_force-van_force)/P_MASS(p);
                }
                else
                        result=0;        // if lift<=adhension, no force added
        }
        else
                result=0; // for cell not neighbor to wall, set to zero       
        return result;
}

83136715

非常感谢你耐心的回复。我想问一下，F_R(wall_f,wall_tf);这个函数是啥意思捏

83136715

我用的是12.1，用你的CODE总是出错捏

83136715

另外你在这个帖子里http://www.cfluid.com/bbs/viewthread.php?tid=120918回复说，F_R（）边界面上是取不到密度值的，可是上面的code为啥用了F_R（）而不是C_R（）呢？