fluent中在动量方程中添加弹性项以及添加5个用户自定义标量函数UDS后
-
我想在fluent中使用黏弹性流体(FENE-P模型),需要使用udf在动量方程中添加弹性项源项,并且使用uds添加6个标量方程(FENE-P模型的本构方程Cxx,Cxy,Cxz,Cyy,Cyz,Czz),与N-S方程耦合求解,编译成功后可以计算,但是uds的值变化很小(初始化时uds的值给的是1,计算结束后基本还是1,知识小数点后几位有变化),几乎体现不出弹性项的作用。求助各位大佬
(1)会不会是多个uds方程需要对i进行定义
在帮助手册找到这个(如果要求解多个标量,则可以在UDF中使用条件IF语句为每个i定义不同的通量函数。i=0与uds-0(要求解的第一个标量方程)相关联。)但是没找到具体怎么操作
FENE-P 模型的本构方程是
直角坐标下为
添加弹性项后的X方向上的动量方程
添加弹性项后的y方向上的动量方程
添加弹性项后的z方向上的动量方程
具体的udf代码#include "mem.h" #define miup (0.00736) #define lamda (2.47024) #define L (40) DEFINE_UDS_UNSTEADY(No_Cxx_unsteady,c,t,i,apu,su) //cxx非稳态项 { real physical_dt, vol, rho, phi_old; physical_dt = RP_Get_Real("physical-time-step"); vol = C_VOLUME(c,t); *apu = -vol / physical_dt;/*implicit part*/ phi_old = C_STORAGE_R(c,t,SV_UDSI_M1(0)); *su = vol*phi_old/physical_dt;/*explicit part*/ } DEFINE_UDS_UNSTEADY(No_Cxy_unsteady,c,t,i,apu,su) //cxy非稳态项 { real physical_dt, vol, rho, phi_old; physical_dt = RP_Get_Real("physical-time-step"); vol = C_VOLUME(c,t); *apu = -vol / physical_dt;/*implicit part*/ phi_old = C_STORAGE_R(c,t,SV_UDSI_M1(1)); *su = vol*phi_old/physical_dt;/*explicit part*/ } DEFINE_UDS_UNSTEADY(No_Cxz_unsteady,c,t,i,apu,su) //cxz非稳态项 { real physical_dt, vol, rho, phi_old; physical_dt = RP_Get_Real("physical-time-step"); vol = C_VOLUME(c,t); *apu = -vol / physical_dt;/*implicit part*/ phi_old = C_STORAGE_R(c,t,SV_UDSI_M1(2)); *su = vol*phi_old/physical_dt;/*explicit part*/ } DEFINE_UDS_UNSTEADY(No_Cyy_unsteady,c,t,i,apu,su) //cyy非稳态项 { real physical_dt, vol, rho, phi_old; physical_dt = RP_Get_Real("physical-time-step"); vol = C_VOLUME(c,t); *apu = -vol / physical_dt;/*implicit part*/ phi_old = C_STORAGE_R(c,t,SV_UDSI_M1(3)); *su = vol*phi_old/physical_dt;/*explicit part*/ } DEFINE_UDS_UNSTEADY(No_Cyz_unsteady,c,t,i,apu,su) //cyz非稳态项 { real physical_dt, vol, rho, phi_old; physical_dt = RP_Get_Real("physical-time-step"); vol = C_VOLUME(c,t); *apu = -vol / physical_dt;/*implicit part*/ phi_old = C_STORAGE_R(c,t,SV_UDSI_M1(4)); *su = vol*phi_old/physical_dt;/*explicit part*/ } DEFINE_UDS_UNSTEADY(No_Czz_unsteady,c,t,i,apu,su) //czz非稳态项 { real physical_dt, vol, rho, phi_old; physical_dt = RP_Get_Real("physical-time-step"); vol = C_VOLUME(c,t); *apu = -vol / physical_dt;/*implicit part*/ phi_old = C_STORAGE_R(c,t,SV_UDSI_M1(5)); *su = vol*phi_old/physical_dt;/*explicit part*/ } DEFINE_UDS_FLUX(No_Cxx_FLUX,f,t,i) //Cxx对流项 { cell_t c0, c1 = -1; Thread *t0, *t1 = NULL; real NV_VEC(psi_vec), NV_VEC(A), flux = 0.0; c0 = F_C0(f,t); t0 = F_C0_THREAD(f,t); F_AREA(A, f, t); if (BOUNDARY_FACE_THREAD_P(t)) { real dens; if (NNULLP(THREAD_STORAGE(t,SV_DENSITY))) dens = F_R(f,t); else dens = C_R(c0,t0); NV_DS(psi_vec, =, F_U(f,t), F_V(f,t), F_W(f,t), *, 1); flux = NV_DOT(psi_vec, A); } else { c1 = F_C1(f,t); t1 = F_C1_THREAD(f,t); NV_DS(psi_vec, =, C_U(c0,t0),C_V(c0,t0),C_W(c0,t0),*,1); NV_DS(psi_vec, +=, C_U(c1,t1),C_V(c1,t1),C_W(c1,t1),*,1); flux = NV_DOT(psi_vec, A)/2.0; } return flux; } DEFINE_UDS_FLUX(No_Cxy_FLUX,f,t,i) //Cxy对流项 { cell_t c0, c1 = -1; Thread *t0, *t1 = NULL; real NV_VEC(psi_vec), NV_VEC(A), flux = 0.0; c0 = F_C0(f,t); t0 = F_C0_THREAD(f,t); F_AREA(A, f, t); if (BOUNDARY_FACE_THREAD_P(t)) { real dens; if (NNULLP(THREAD_STORAGE(t,SV_DENSITY))) dens = F_R(f,t); else dens = C_R(c0,t0); NV_DS(psi_vec, =, F_U(f,t), F_V(f,t), F_W(f,t), *, 1); flux = NV_DOT(psi_vec, A); } else { c1 = F_C1(f,t); t1 = F_C1_THREAD(f,t); NV_DS(psi_vec, =, C_U(c0,t0),C_V(c0,t0),C_W(c0,t0),*,1); NV_DS(psi_vec, +=, C_U(c1,t1),C_V(c1,t1),C_W(c1,t1),*,1); flux = NV_DOT(psi_vec, A)/2.0; } return flux; } DEFINE_UDS_FLUX(No_Cxz_FLUX,f,t,i) //Cxz对流项 { cell_t c0, c1 = -1; Thread *t0, *t1 = NULL; real NV_VEC(psi_vec), NV_VEC(A), flux = 0.0; c0 = F_C0(f,t); t0 = F_C0_THREAD(f,t); F_AREA(A, f, t); if (BOUNDARY_FACE_THREAD_P(t)) { real dens; if (NNULLP(THREAD_STORAGE(t,SV_DENSITY))) dens = F_R(f,t); else dens = C_R(c0,t0); NV_DS(psi_vec, =, F_U(f,t), F_V(f,t), F_W(f,t), *, 1); flux = NV_DOT(psi_vec, A); } else { c1 = F_C1(f,t); t1 = F_C1_THREAD(f,t); NV_DS(psi_vec, =, C_U(c0,t0),C_V(c0,t0),C_W(c0,t0),*,1); NV_DS(psi_vec, +=, C_U(c1,t1),C_V(c1,t1),C_W(c1,t1),*,1); flux = NV_DOT(psi_vec, A)/2.0; } return flux; } DEFINE_UDS_FLUX(No_Cyy_FLUX,f,t,i) //Cyy对流项 { cell_t c0, c1 = -1; Thread *t0, *t1 = NULL; real NV_VEC(psi_vec), NV_VEC(A), flux = 0.0; c0 = F_C0(f,t); t0 = F_C0_THREAD(f,t); F_AREA(A, f, t); if (BOUNDARY_FACE_THREAD_P(t)) { real dens; if (NNULLP(THREAD_STORAGE(t,SV_DENSITY))) dens = F_R(f,t); else dens = C_R(c0,t0); NV_DS(psi_vec, =, F_U(f,t), F_V(f,t), F_W(f,t), *, 1); flux = NV_DOT(psi_vec, A); } else { c1 = F_C1(f,t); t1 = F_C1_THREAD(f,t); NV_DS(psi_vec, =, C_U(c0,t0),C_V(c0,t0),C_W(c0,t0),*,1); NV_DS(psi_vec, +=, C_U(c1,t1),C_V(c1,t1),C_W(c1,t1),*,1); flux = NV_DOT(psi_vec, A)/2.0; } return flux; } DEFINE_UDS_FLUX(No_Cyz_FLUX,f,t,i) //Cyz对流项 { cell_t c0, c1 = -1; Thread *t0, *t1 = NULL; real NV_VEC(psi_vec), NV_VEC(A), flux = 0.0; c0 = F_C0(f,t); t0 = F_C0_THREAD(f,t); F_AREA(A, f, t); if (BOUNDARY_FACE_THREAD_P(t)) { real dens; if (NNULLP(THREAD_STORAGE(t,SV_DENSITY))) dens = F_R(f,t); else dens = C_R(c0,t0); NV_DS(psi_vec, =, F_U(f,t), F_V(f,t), F_W(f,t), *, 1); flux = NV_DOT(psi_vec, A); } else { c1 = F_C1(f,t); t1 = F_C1_THREAD(f,t); NV_DS(psi_vec, =, C_U(c0,t0),C_V(c0,t0),C_W(c0,t0),*,1); NV_DS(psi_vec, +=, C_U(c1,t1),C_V(c1,t1),C_W(c1,t1),*,1); flux = NV_DOT(psi_vec, A)/2.0; } return flux; } DEFINE_UDS_FLUX(No_Czz_FLUX,f,t,i) //Czz对流项 { cell_t c0, c1 = -1; Thread *t0, *t1 = NULL; real NV_VEC(psi_vec), NV_VEC(A), flux = 0.0; c0 = F_C0(f,t); t0 = F_C0_THREAD(f,t); F_AREA(A, f, t); if (BOUNDARY_FACE_THREAD_P(t)) { real dens; if (NNULLP(THREAD_STORAGE(t,SV_DENSITY))) dens = F_R(f,t); else dens = C_R(c0,t0); NV_DS(psi_vec, =, F_U(f,t), F_V(f,t), F_W(f,t), *, 1); flux = NV_DOT(psi_vec, A); } else { c1 = F_C1(f,t); t1 = F_C1_THREAD(f,t); NV_DS(psi_vec, =, C_U(c0,t0),C_V(c0,t0),C_W(c0,t0),*,1); NV_DS(psi_vec, +=, C_U(c1,t1),C_V(c1,t1),C_W(c1,t1),*,1); flux = NV_DOT(psi_vec, A)/2.0; } return flux; } DEFINE_SOURCE(No_Cxx_source, c, t, dS, eqn) //Cxx源项 { real source; real fR; fR=(L*L-3)/(L*L-C_UDSI(c,t,0)-C_UDSI(c,t,3)-C_UDSI(c,t,5)); source=2*(C_UDSI(c,t,0)*C_DUDX(c,t)+C_UDSI(c,t,1)*C_DUDY(c,t)+C_UDSI(c,t,2)*C_DUDZ(c,t))-(fR*C_UDSI(c,t,0)-1)/lamda; dS[eqn]=0.0; return source; } DEFINE_SOURCE(No_Cxy_source, c, t, dS, eqn) //Cxy源项 { real source; real fR; fR=(L*L-3)/(L*L-C_UDSI(c,t,0)-C_UDSI(c,t,3)-C_UDSI(c,t,5)); source=(C_UDSI(c,t,0)*C_DVDX(c,t)+C_UDSI(c,t,1)*C_DVDY(c,t)+C_UDSI(c,t,2)*C_DVDZ(c,t)+\ C_UDSI(c,t,1)*C_DUDX(c,t)+C_UDSI(c,t,3)*C_DUDY(c,t)+C_UDSI(c,t,4)*C_DUDZ(c,t))-\ (fR*C_UDSI(c,t,1))/lamda; dS[eqn]=0.0; return source; } DEFINE_SOURCE(No_Cxz_source, c, t, dS, eqn) //Cxz源项 { real source; real fR; fR=(L*L-3)/(L*L-C_UDSI(c,t,0)-C_UDSI(c,t,3)-C_UDSI(c,t,5)); source=(C_UDSI(c,t,0)*C_DWDX(c,t)+C_UDSI(c,t,1)*C_DWDY(c,t)+C_UDSI(c,t,2)*C_DWDZ(c,t)+\ C_UDSI(c,t,2)*C_DUDX(c,t)+C_UDSI(c,t,4)*C_DUDY(c,t)+C_UDSI(c,t,5)*C_DUDZ(c,t))-\ (fR*C_UDSI(c,t,2))/lamda; dS[eqn]=0.0; return source; } DEFINE_SOURCE(No_Cyy_source, c, t, dS, eqn) //Cyy源项 { real source; real fR; fR=(L*L-3)/(L*L-C_UDSI(c,t,0)-C_UDSI(c,t,3)-C_UDSI(c,t,5)); source=2*(C_UDSI(c,t,1)*C_DVDX(c,t)+C_UDSI(c,t,3)*C_DVDY(c,t)+C_UDSI(c,t,4)*C_DVDZ(c,t))-(fR*C_UDSI(c,t,3)-1)/lamda; dS[eqn]=0.0; return source; } DEFINE_SOURCE(No_Cyz_source, c, t, dS, eqn) //Cyz源项 { real source; real fR; fR=(L*L-3)/(L*L-C_UDSI(c,t,0)-C_UDSI(c,t,3)-C_UDSI(c,t,5)); source=(C_UDSI(c,t,1)*C_DWDX(c,t)+C_UDSI(c,t,3)*C_DWDY(c,t)+C_UDSI(c,t,4)*C_DWDZ(c,t)+\ C_UDSI(c,t,2)*C_DVDX(c,t)+C_UDSI(c,t,4)*C_DVDY(c,t)+C_UDSI(c,t,5)*C_DVDZ(c,t))-\ (fR*C_UDSI(c,t,4))/lamda; dS[eqn]=0.0; return source; } DEFINE_SOURCE(No_Czz_source, c, t, dS, eqn) //Czz源项 { real source; real fR; fR=(L*L-3)/(L*L-C_UDSI(c,t,0)-C_UDSI(c,t,3)-C_UDSI(c,t,5)); source=2*(C_UDSI(c,t,2)*C_DWDX(c,t)+C_UDSI(c,t,4)*C_DWDY(c,t)+C_UDSI(c,t,5)*C_DWDZ(c,t))-(fR*C_UDSI(c,t,5)-1)/lamda; dS[eqn]=0.0; return source; } DEFINE_SOURCE(x_momentum, c, t, dS, eqn) //x方向动量源项 { real source; real cxx = C_UDSI(c,t,0); real cxy = C_UDSI(c,t,1); real cxz = C_UDSI(c,t,2); real cyx = C_UDSI(c,t,1); real cyy = C_UDSI(c,t,3); real cyz = C_UDSI(c,t,4); real czx = C_UDSI(c,t,2); real czy = C_UDSI(c,t,4); real czz = C_UDSI(c,t,5); real NV_VEC(gcxx); real NV_VEC(gcxy); real NV_VEC(gcxz); real NV_VEC(gcyx); real NV_VEC(gcyy); real NV_VEC(gcyz); real NV_VEC(gczx); real NV_VEC(gczy); real NV_VEC(gczz); NV_DS(gcxx, =,C_UDSI_G(c,t,0)[0],C_UDSI_G(c,t,0)[1],C_UDSI_G(c,t,0)[2],*,1); NV_DS(gcxy, =,C_UDSI_G(c,t,1)[0],C_UDSI_G(c,t,1)[1],C_UDSI_G(c,t,1)[2],*,1); NV_DS(gcxz, =,C_UDSI_G(c,t,2)[0],C_UDSI_G(c,t,2)[1],C_UDSI_G(c,t,2)[2],*,1); NV_DS(gcyx, =,C_UDSI_G(c,t,1)[0],C_UDSI_G(c,t,1)[1],C_UDSI_G(c,t,1)[2],*,1); NV_DS(gcyy, =,C_UDSI_G(c,t,3)[0],C_UDSI_G(c,t,3)[1],C_UDSI_G(c,t,3)[2],*,1); NV_DS(gcyz, =,C_UDSI_G(c,t,4)[0],C_UDSI_G(c,t,4)[1],C_UDSI_G(c,t,4)[2],*,1); NV_DS(gczx, =,C_UDSI_G(c,t,2)[0],C_UDSI_G(c,t,2)[1],C_UDSI_G(c,t,2)[2],*,1); NV_DS(gczy, =,C_UDSI_G(c,t,4)[0],C_UDSI_G(c,t,4)[1],C_UDSI_G(c,t,4)[2],*,1); NV_DS(gczz, =,C_UDSI_G(c,t,5)[0],C_UDSI_G(c,t,5)[1],C_UDSI_G(c,t,5)[2],*,1); source = miup*(L*L-3)/lamda/pow(L*L-cxx-cyy-czz,2)*((gcxx[0]+gcxy[1]+gcxz[2])*(L*L-cxx-cyy-czz)+\ cxx*(gcxx[0]+gcyy[0]+gczz[0])+cxy*(gcxx[1]+gcyy[1]+gczz[1])+cxz*(gcxx[2]+gcyy[2]+gczz[2])); return source; } DEFINE_SOURCE(y_momentum, c, t, dS, eqn) //y方向动量源项 { real source; real cxx = C_UDSI(c,t,0); real cxy = C_UDSI(c,t,1); real cxz = C_UDSI(c,t,2); real cyx = C_UDSI(c,t,1); real cyy = C_UDSI(c,t,3); real cyz = C_UDSI(c,t,4); real czx = C_UDSI(c,t,2); real czy = C_UDSI(c,t,4); real czz = C_UDSI(c,t,5); real NV_VEC(gcxx); real NV_VEC(gcxy); real NV_VEC(gcxz); real NV_VEC(gcyx); real NV_VEC(gcyy); real NV_VEC(gcyz); real NV_VEC(gczx); real NV_VEC(gczy); real NV_VEC(gczz); NV_DS(gcxx, =,C_UDSI_G(c,t,0)[0],C_UDSI_G(c,t,0)[1],C_UDSI_G(c,t,0)[2],*,1); NV_DS(gcxy, =,C_UDSI_G(c,t,1)[0],C_UDSI_G(c,t,1)[1],C_UDSI_G(c,t,1)[2],*,1); NV_DS(gcxz, =,C_UDSI_G(c,t,2)[0],C_UDSI_G(c,t,2)[1],C_UDSI_G(c,t,2)[2],*,1); NV_DS(gcyx, =,C_UDSI_G(c,t,1)[0],C_UDSI_G(c,t,1)[1],C_UDSI_G(c,t,1)[2],*,1); NV_DS(gcyy, =,C_UDSI_G(c,t,3)[0],C_UDSI_G(c,t,3)[1],C_UDSI_G(c,t,3)[2],*,1); NV_DS(gcyz, =,C_UDSI_G(c,t,4)[0],C_UDSI_G(c,t,4)[1],C_UDSI_G(c,t,4)[2],*,1); NV_DS(gczx, =,C_UDSI_G(c,t,2)[0],C_UDSI_G(c,t,2)[1],C_UDSI_G(c,t,2)[2],*,1); NV_DS(gczy, =,C_UDSI_G(c,t,4)[0],C_UDSI_G(c,t,4)[1],C_UDSI_G(c,t,4)[2],*,1); NV_DS(gczz, =,C_UDSI_G(c,t,5)[0],C_UDSI_G(c,t,5)[1],C_UDSI_G(c,t,5)[2],*,1); source = miup*(L*L-3)/lamda/pow(L*L-cxx-cyy-czz,2)*((gcyx[0]+gcyy[1]+gcyz[2])*(L*L-cxx-cyy-czz)+\ cyx*(gcxx[0]+gcyy[0]+gczz[0])+cyy*(gcxx[1]+gcyy[1]+gczz[1])+cyz*(gcxx[2]+gcyy[2]+gczz[2])); return source; } DEFINE_SOURCE(z_momentum, c, t, dS, eqn) //z方向动量源项 { real source; real cxx = C_UDSI(c,t,0); real cxy = C_UDSI(c,t,1); real cxz = C_UDSI(c,t,2); real cyx = C_UDSI(c,t,1); real cyy = C_UDSI(c,t,3); real cyz = C_UDSI(c,t,4); real czx = C_UDSI(c,t,2); real czy = C_UDSI(c,t,4); real czz = C_UDSI(c,t,5); real NV_VEC(gcxx); real NV_VEC(gcxy); real NV_VEC(gcxz); real NV_VEC(gcyx); real NV_VEC(gcyy); real NV_VEC(gcyz); real NV_VEC(gczx); real NV_VEC(gczy); real NV_VEC(gczz); NV_DS(gcxx, =,C_UDSI_G(c,t,0)[0],C_UDSI_G(c,t,0)[1],C_UDSI_G(c,t,0)[2],*,1); NV_DS(gcxy, =,C_UDSI_G(c,t,1)[0],C_UDSI_G(c,t,1)[1],C_UDSI_G(c,t,1)[2],*,1); NV_DS(gcxz, =,C_UDSI_G(c,t,2)[0],C_UDSI_G(c,t,2)[1],C_UDSI_G(c,t,2)[2],*,1); NV_DS(gcyx, =,C_UDSI_G(c,t,1)[0],C_UDSI_G(c,t,1)[1],C_UDSI_G(c,t,1)[2],*,1); NV_DS(gcyy, =,C_UDSI_G(c,t,3)[0],C_UDSI_G(c,t,3)[1],C_UDSI_G(c,t,3)[2],*,1); NV_DS(gcyz, =,C_UDSI_G(c,t,4)[0],C_UDSI_G(c,t,4)[1],C_UDSI_G(c,t,4)[2],*,1); NV_DS(gczx, =,C_UDSI_G(c,t,2)[0],C_UDSI_G(c,t,2)[1],C_UDSI_G(c,t,2)[2],*,1); NV_DS(gczy, =,C_UDSI_G(c,t,4)[0],C_UDSI_G(c,t,4)[1],C_UDSI_G(c,t,4)[2],*,1); NV_DS(gczz, =,C_UDSI_G(c,t,5)[0],C_UDSI_G(c,t,5)[1],C_UDSI_G(c,t,5)[2],*,1); source = miup*(L*L-3)/lamda/pow(L*L-cxx-cyy-czz,2)*((gczx[0]+gczy[1]+gczz[2])*(L*L-cxx-cyy-czz)+\ czx*(gcxx[0]+gcyy[0]+gczz[0])+czy*(gcxx[1]+gcyy[1]+gczz[1])+czz*(gcxx[2]+gcyy[2]+gczz[2])); return source; } code_text
-
@bestucan 谢谢老师,我还想问一下我uds方程里对流项没有密度,我直接把udf帮助文档里面对流项的代码中密度的地方改为1可以吧?
DEFINE_UDS_FLUX(No_Cxx_FLUX,f,t,i) //Cxx对流项 { cell_t c0, c1 = -1; Thread *t0, *t1 = NULL; real NV_VEC(psi_vec), NV_VEC(A), flux = 0.0; c0 = F_C0(f,t); t0 = F_C0_THREAD(f,t); F_AREA(A, f, t); if (BOUNDARY_FACE_THREAD_P(t)) { real dens; if (NNULLP(THREAD_STORAGE(t,SV_DENSITY))) dens = F_R(f,t); else dens = C_R(c0,t0); NV_DS(psi_vec, =, F_U(f,t), F_V(f,t), F_W(f,t), *, 1); flux = NV_DOT(psi_vec, A); } else { c1 = F_C1(f,t); t1 = F_C1_THREAD(f,t); NV_DS(psi_vec, =, C_U(c0,t0),C_V(c0,t0),C_W(c0,t0),*,1); NV_DS(psi_vec, +=, C_U(c1,t1),C_V(c1,t1),C_W(c1,t1),*,1); flux = NV_DOT(psi_vec, A)/2.0; } return flux; }