求助：DEFINE_ADJUST求两相界面力能load但不能计算

guohf00001

我用vof法计算气液两相流，想先找出两相界面，并求界面电密度，用DEFINE_ADJUST和DEFINE_ON_DEMAND都试过，能顺利load进去，但一hook这个界面电密度文件fluent就提示致命错误。下面是我的部分程序，请高手指导一下，先谢了。

为方便，附上附件，最后一个udf DEFINE_ADJUST(calm_rhos, domain)出错。   vof_interface_force.doc (3.23 KB, 下载次数: 3)

2017-2-26 12:02 上传

点击文件名下载附件

#include "stdio.h"
#include "udf.h"
#include "sg.h"
#include "sg_mphase.h"
#include "math.h"
#include "flow.h"
#include "mem.h"
#include "metric.h"
#include "id.h"

#define EPS_1 4.25   /* 相对介电常数*/
#define EPS_0 8.8541878e-12  /* electrical permittivity 真空介电常数c^2/Nm^2*/


/* Define which user-defined scalars to use. */

enum                 /*枚举类型名*/
{
        phi,    /* 电势*/  

        N_REQUIRED_UDSI    /*存取用户自定义标量输运方程的数目*/
};


enum                 /*枚举类型名*/
{
        Ex,
        Ey,
        Ez,
        vfx,
        vfy,
        vfz,
        rhos,  /* 界面电密度*/
        N_REQUIRED_UDM    /*存取用户自定义标量输运方程的数目*/
};

DEFINE_ADJUST(store_VOF_gradient, d)  /*存储vof体积分数的梯度*/
{
        Thread *t;
        Thread *ppt;
        Thread **pt;
        cell_t c;
        int phase_domain_index=1;
        Domain *pDomain = DOMAIN_SUB_DOMAIN(d,phase_domain_index);
        Alloc_Storage_Vars(pDomain,SV_VOF_RG,SV_VOF_G,SV_NULL);
        Scalar_Reconstruction(pDomain, SV_VOF,-1,SV_VOF_RG,NULL);
        Scalar_Derivatives(pDomain,SV_VOF,-1,SV_VOF_G,SV_VOF_RG, Vof_Deriv_Accumulate);

        mp_thread_loop_c(t,d,pt)
        {
                if (FLUID_THREAD_P(t))
                {
                        ppt = pt[phase_domain_index];

                        begin_c_loop (c,t)
                        {
                                C_UDMI(c,t,vfx) = C_VOF_G(c,ppt)[0];
                                C_UDMI(c,t,vfy) = C_VOF_G(c,ppt)[1];
                                C_UDMI(c,t,vfz) = C_VOF_G(c,ppt)[2];
                        }
                        end_c_loop (c,t)
                }
        }

        Free_Storage_Vars(pDomain,SV_VOF_RG,SV_VOF_G,SV_NULL);
}


DEFINE_ADJUST(Calc_E_xyz, domain)  /*存储电场强度*/
{

        Thread *t;
        cell_t c;
        face_t f;
        Thread *t0;
        cell_t c0;
        real dr0[ND_ND],dr1[ND_ND], xf[ND_ND], dy;


        if (n_uds < N_REQUIRED_UDSI)

                Internal_Error("not enough user-defined scalars allocated");


        if(n_udm < N_REQUIRED_UDM)

                Error("Not enough udm allocated\n");


        /* Do nothing if gradient isn't allocated yet. */

        if (!Data_Valid_P()) return;

        thread_loop_c(t,domain)
        {
                if (FLUID_THREAD_P(t))
                        if (NULL != THREAD_STORAGE(t,SV_UDS_I(phi)) && NULL != T_STORAGE_R_NV (t,SV_UDSI_G(phi)))
                        {
                                begin_c_loop(c,t)       /*体循环*/
                                {
                                        C_UDMI(c,t,Ex) = -1.*C_UDSI_G(c,t,phi)[0];    /*Ex*/
                                        C_UDMI(c,t,Ey) = -1.*C_UDSI_G(c,t,phi)[1];
                                        C_UDMI(c,t,Ez) = -1.*C_UDSI_G(c,t,phi)[2];
                                }
                                end_c_loop(c,t)
                                        /* Message("x electric strength: %g\n", C_UDMI(c,t,Ex)); */
                        }
        }

}


DEFINE_ADJUST(calm_rhos, domain)   /* 界面电密度，出问题的地方*/
{       
        Thread *t;
        Thread *pri_th,*sec_th;
        cell_t c;
        real xc[ND_ND], EEp[3], EEs[3], nn[3], EE[3], rr;
        /*  Domain *domain = Get_Domain(1); Get domain pointer */
        int ID = 1;
        t =Lookup_Thread(domain, ID); /* mixture-level thread pointer */

        pri_th = THREAD_SUB_THREAD(t,0);
        sec_th = THREAD_SUB_THREAD(t,1);
       
        thread_loop_c(t,domain)
        {
                begin_c_loop(c,t)       /*体循环*/
                {
                        C_CENTROID(xc,c,t);
                        if (C_VOF(c,sec_th)>0. && C_VOF(c,sec_th)<1.)
                        {
                                EEs[0]=C_UDMI(c,sec_th,Ex);/* 次相电场强度 */
                                EEs[1]=C_UDMI(c,sec_th,Ey);
                                EEs[2]=C_UDMI(c,sec_th,Ez);
                                EEp[0]=C_UDMI(c,pri_th,Ex);/* 基相电场强度 */
                                EEp[1]=C_UDMI(c,pri_th,Ey);       
                                EEp[2]=C_UDMI(c,pri_th,Ez);
                                nn[0]=C_UDMI(c,sec_th,vfx);/* vof梯度、法向矢量 */
                                nn[1]=C_UDMI(c,sec_th,vfy);
                                nn[2]=C_UDMI(c,sec_th,vfz);
                                NV_V(EEs,=,EEp);
                                rr=EPS_0*EPS_1;
                                NV_VS_VS(EE, =,EEp,*,EPS_0, -, EEs,*,rr);
                                C_UDMI(c,t,rhos) = NV_DOT(nn,EE);
                        }
                        else
                        {
                                C_UDMI(c,t,rhos) =0;
                        }
                }
                end_c_loop(c,t)
        }
}

guohf00001

没人理啊，急啊

aaa-1234

可以把case一起发给我，我帮你看看，挺不错的东西
加我QQ 103614652

guohf00001

aaa-1234 发表于 2017-3-1 00:55
可以把case一起发给我，我帮你看看，挺不错的东西
加我QQ 103614652

非常感谢

guohf00001

aaa-1234 发表于 2017-3-1 00:55
可以把case一起发给我，我帮你看看，挺不错的东西
加我QQ 103614652

我的QQ 247136248。我已加你为好友了，谢谢帮忙。