我想在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 模型的本构方程是
1db2aecf-bb28-4e18-83b2-b78012d300a0-8de57d7561f6f0ce3fa3533604fefd3.png
直角坐标下为
2c357f13-0b20-4832-8817-2d63da569bd0-4ad7186b8fd9973a10024ee26d56f59.png
49efe7aa-242a-40a1-99f0-77492f023026-e3f29d15bfaa922b3ad458e6135be3c.png
2f7eef5f-0dc2-481c-ae2e-fbf59aaf90f4-2675ed868f2111033ffc1d7af482a07.png
添加弹性项后的X方向上的动量方程
7da0ce9c-10e5-467d-bbc7-a5e298f58678-8985c92d1edb7e6c3b4912f0ccb09ae.png
添加弹性项后的y方向上的动量方程
3c27a4d7-a3de-4795-97bf-85f89b1ca171-7f1fc03c84026987030e60dd8546d30.png
添加弹性项后的z方向上的动量方程
3ba2dfa4-8db9-4ae7-ac16-e72c39f52553-4506fbbd7bf7138e779729d1836a718.png
具体的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
CFD中文网

冠竹
