CFD中文网

    CFD中文网

    • 登录
    • 搜索
    • 最新

    UDS方程计算出来的污染物总量不守恒

    Fluent
    3
    11
    1582
    正在加载更多帖子
    • 从旧到新
    • 从新到旧
    • 最多赞同
    回复
    • 在新帖中回复
    登录后回复
    此主题已被删除。只有拥有主题管理权限的用户可以查看。
    • 深蓝
      深蓝 最后由 李东岳 编辑

      问题: UDS方程在封闭计算域内的计算出来的污染物总量不守恒。
      详细描述:
      我需要用FLUENT的UDS方程计算污染物的传播
      UDS标量方法:
      $$\frac { \partial \rho \phi } { \partial t } + \frac { \partial \rho U \phi} { \partial x } = \frac { \partial } { \partial x } (\Gamma _ { \Phi ,eff } \frac { \partial \phi } { \partial x } ) + S _ { 0 }$$
      我在一个正方形内计算污染物的扩散,无进出口,长宽为10m* 10m,具体如图:
      图片2.png
      在计算域中间设置了UDS的源项,源的大小为0.2m*0.2m,释放时间为0.1s,代码很简单,如下:
      第一种方法:

      DEFINE_SOURCE(aerosol_source,c,t,dS,eqn)
      {
         real x[ND_ND];
         real con, source;
         C_CENTROID(x,c,t);
        
         if(x[0]>4.5&&x[0]<5.0&&x[1]>4.5&&x[1]<5.0)
             {
              source=1;
              }
         else
             {
              source=0;
              }
      
         dS[eqn] = 0;
         return source;
      }
      
      

      第二种方法
      或者Mark一个区域出来,把这个区域的source设为1,其他区域的source设为0。如图:
      图片8.png

      所说的守恒就是源的释放量(udf手册说DEFINE_SOURCE宏的source的单位是产生率/m3,如果是质量源就是kg/(s-m3))乘释放时间,再乘释放体积应该等于总量,但是上述两种source设置方法得到的结果在计算域中统计都不守恒。
      这两种方法计算出来结果还不同,如下:
      第一种方法:
      图片3.png
      图片5.png
      第二种方法:
      图片4.png
      图片6.png
      第一种方法为15.99,第二种方法为2.04。
      疑问:
      按理来说应该等于$0.2* 0.2*0.1=0.004$,这算出来跟释放源完全不守恒,请问大家这是什么原因呢?:135:

      1 条回复 最后回复 回复 引用
      • bestucan
        bestucan 版主 副教授 最后由 编辑

        那个设置的值就是“值”,不是“通量”。就好比,设定一个一百度的铁球在那里,但是不是说把一百度铁球上的热量全都放到场内。

        如果环境温度800度,这个100度的铁球还会吸收热量。吸收了应该升温,但是每次迭代都会又设定为100度。所以不停地吸收。

        如果环境温度0度,这个100度的铁球持续释放热量。释放了应该降温,但是每次迭代都会又设定为100度,所以不停的释放。

        是否吸收释放,以及吸收释放的速度取决于温度差。这里应该是浓度差。时间步越小,越能一直保持较高的浓度差,释放的就越多

        就好比,给了个压力入口边界条件,却想得到固定的通量。

        滚来滚去……~(~o ̄▽ ̄)~o 滚来滚去都不能让大家看出来我不是老师么 O_o

        异步沟通方式(《posting style》from wiki)(下载后打开):
        https://www.jianguoyun.com/p/Dc52X2sQsLv2BRiqnKYD
        提问的智慧(github在gitee的镜像):
        https://gitee.com/bestucan/How-To-Ask-Questions-The-Smart-Way/blob/master/README-zh_CN.md

        深蓝 1 条回复 最后回复 回复 引用
        • 深蓝
          深蓝 @bestucan 最后由 编辑

          @bestucan 谢谢您得回复,在fluent 的udf手册中写道,源项的单位是generation-rate/volume,单位体积的生成率在往外传递的时候不应该是一个通量吗?设置的“值”或是设置“通量”是否是对应的边界条件,而这个地方设置的是源项。

          如果我的目的是放置一个1kg的污染物在这儿,看他随时间扩散后质量的变化情况,请问改怎么设置啊。
          818.png

          bestucan 1 条回复 最后回复 回复 引用
          • bestucan
            bestucan 版主 副教授 @深蓝 最后由 编辑

            @深蓝 可以做这样一个测试,全部区域都给这个源,看浓度是不是线性升高的就知道了是哪种了。

            如果只是放置个污染源看扩散,就在计算初始化的时候设定。这样应该比较简单。
            大部分计算初始化都是allzone 设定为 0 。可以在patch里分区域设定计算初始状态。
            面积乘以UDS值,就是初始污染物的量。想要单位1的污染物。就给单位面积上设定uds的初始值为1。

            滚来滚去……~(~o ̄▽ ̄)~o 滚来滚去都不能让大家看出来我不是老师么 O_o

            异步沟通方式(《posting style》from wiki)(下载后打开):
            https://www.jianguoyun.com/p/Dc52X2sQsLv2BRiqnKYD
            提问的智慧(github在gitee的镜像):
            https://gitee.com/bestucan/How-To-Ask-Questions-The-Smart-Way/blob/master/README-zh_CN.md

            深蓝 1 条回复 最后回复 回复 引用
            • 深蓝
              深蓝 @bestucan 最后由 编辑

              @bestucan 谢谢您得回复,通过您得建议,经过我几次尝试,结果如下:

              方法一:
              初始化所有区域的uds=0,u=2m/s,空气密度为1kg/m3
              Patch方块的uds=1
              方块的网格数量为29*29
              uds的总和为 29*29*1=841
              经过0.5s的发展如右图,uds的总和还是841
              Note:
              1、如果密度为2 kg/m3,则uds的总和为420
              8191.png

              方法二:
              初始化所有区域的uds=0,空气密度为1kg/m3
              设置方块的source=1
              方块的网格数量为29*29
              进过0.5s的释放和发展如右图(线性增加),uds的总和为420
              Note:
              1、 如果密度为0.5 kg/m3,则uds的总和为840
              2、 如果是稳态计算,uds是线性增加的
              8192.png
              8193.png
              所以,source的单位应该是generation-rate*density/cell numbers

              这里面需要计算网格的数量,我想udf怎么得到网格的数量呢,找了半天没有合适的,C_NNODES可以用来计算cell的节点数量,F_NNODES可以计算网格面上节点的数量,C_NFACES可以计算网格中面的数量,是否能间接得到网格的数量,但是udf手册中没有这个实例,但是在计算的时候编译一直通过不了,呜呜呜:135: ,代码如下:

              #include "udf.h"
              #include "mem.h"
              DEFINE_ON_DEMAND(name)
              {	
              	int sum=0;
              	Thread *t;
              	cell_t c;	
              		begin_c_loop(c, t)
              			{									
              			  {
              			    sum+=C_NNODES(c, t);
              			  }
              			}
              		end_c_loop(c, t)
              			
              	Message0("\n num=%f\n",num/4.0);
              		
              }
              
              
              bestucan 1 条回复 最后回复 回复 引用
              • bestucan
                bestucan 版主 副教授 @深蓝 最后由 编辑

                @深蓝 这个网格数量……我猜应该是网格面积。这里刚好你的网格都一样大,不信你换个网格密度试试:chouchou:

                这个材料密度,按说uds就是在场内运输个标量,跟其他的不太相关,但是一般应用中都会是某种物质,和其他物质有相互作用,然后就得有质量属性。这里材料密度能影响到uds的总量,应该是有个默认的1:1的关系,定义uds那里有个flux function,应该选的是mass flow rate。感兴趣可以找例子在udf写个flux function,兴许uds总量和材料密度的关系就没了。

                这patch都解决问题了,怎么还倒腾udf:qichuang:

                滚来滚去……~(~o ̄▽ ̄)~o 滚来滚去都不能让大家看出来我不是老师么 O_o

                异步沟通方式(《posting style》from wiki)(下载后打开):
                https://www.jianguoyun.com/p/Dc52X2sQsLv2BRiqnKYD
                提问的智慧(github在gitee的镜像):
                https://gitee.com/bestucan/How-To-Ask-Questions-The-Smart-Way/blob/master/README-zh_CN.md

                深蓝 1 条回复 最后回复 回复 引用
                • 深蓝
                  深蓝 @bestucan 最后由 深蓝 编辑

                  @bestucan 谢谢您得持续回复。

                  我进行了三次尝试,200*200=40000个网格,密度=1.0时,整个计算域都添加源项,1.0s后,总量为40000,如果是300*300=90000个网格,其他不变时,总量为90000,如果密度变成2.0,总量就是45000。可见跟密度有有关系的,在UDS对流项的宏(DEFINE_UDS_FLUX)里面需要获取网格的密度,UDF手册的表述如下:
                  未命名1630154273.png

                  我把uds对流项中的密度都设置成1.0常数,没有效果。
                  我尝试把uds扩散项设置成常数也没有效果,分析应该是DEFINE_SOURCE宏返回值是一个扩散系数,密度项已经包含在默认的计算中了。
                  对流项代码如下:

                  DEFINE_UDS_FLUX(QR_flux,f,t,i)
                  {
                      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 face lies at domain boundary, use face values. */
                  /* BOUNDARY_FACE_THREAD_P(t) expands to a function that returns TRUE if Thread *t is a boundary face thread. */
                  
                    	if(BOUNDARY_FACE_THREAD_P(t))
                      	{
                        		real dens;
                  	  
                  /* Depending on its BC, density may not be set on face thread. */
                  /* Test whether the memory for SV_DENSITY has already been allocated on a given Thread or not. */
                  
                        		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), *, dens);
                  
                        		flux = 1 * NV_DOT(psi_vec, A);
                      	}
                    
                  /* If face lies IN the domain, use average of adjacent cells. */
                    
                    	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),*,C_R(c0,t0));
                        		NV_DS(psi_vec, +=, C_U(c1,t1),C_V(c1,t1),C_W(c1,t1),*,C_R(c1,t1));
                  
                        		flux = 1 * NV_DOT(psi_vec, A)/2.0;
                      	}
                  
                    	return flux;
                      
                  }
                  

                  ###另外还有个小问题想请教您一下:
                  未命名1630154958.png
                  问题:UDS的方程中,瞬态项和对流项的处理都有现成可用的宏,我的方程中需要添加一个四阶导数项,看来只能放在源项中,源项处理四阶导数项的话,需要自己编写四阶导数吗,源项的线性化怎么处理……,您对处理这个四阶导数项有什么建议:chouchou:

                  再次谢谢您得回复!!:xinxin2:
                  尝试的代码如下:

                  
                  DEFINE_SOURCE(QR_source_2d, c, t, dS, eqn)
                  {
                  	real source;
                  	real epsilon; 														/*stability coefficient*/
                  	real rho;
                  
                  		
                  C_UDSI(c,t,1) = C_UDSI_G(c,t,0)[0];							/*DcDx, returns the x-component of the UDS gradient vector*/	
                  C_UDSI(c,t,2) = C_UDSI_G(c,t,1)[0];							/*DcDxx, returns the x-component of the UDS gradient vector*/	
                  C_UDSI(c,t,3) = C_UDSI_G(c,t,2)[0];							/*DcDxxx, returns the x-component of the UDS gradient vector*/	
                  C_UDSI(c,t,4) = C_UDSI_G(c,t,3)[0];							/*DcDxxxx, returns the x-component of the UDS gradient vector*/		
                  C_UDSI(c,t,5) = C_UDSI_G(c,t,0)[1];							/*DcDy, returns the y-component of the UDS gradient vector*/	
                  C_UDSI(c,t,6) = C_UDSI_G(c,t,5)[1];							/*DcDy, returns the y-component of the UDS gradient vector*/	
                  C_UDSI(c,t,7) = C_UDSI_G(c,t,6)[1];							/*DcDyyy, returns the y-component of the UDS gradient vector*/	
                  C_UDSI(c,t,8) = C_UDSI_G(c,t,7)[1];							/*DcDyyyy, returns the y-component of the UDS gradient vector*/
                  
                  epsilon = 5*pow(10,-6);
                  rho = C_R(c,t);
                  source = rho * epsilon * (C_UDSI(c,t,4) + C_UDSI(c,t,8));	/*dcdxxxx+dcdyyyy*/
                  			
                  dS[eqn] = 0;												/*source linearization*/
                  			
                  return source;
                  					
                  }
                  
                  
                  
                  李东岳 1 条回复 最后回复 回复 引用
                  • 李东岳
                    李东岳 管理员 @深蓝 最后由 编辑

                    @深蓝 插一句,四阶导数项应该就是对$k=\nabla\cdot\nabla \phi$再做拉普拉斯:$\nabla\cdot\nabla k$

                    线上CFD课程开始报名:http://www.dyfluid.com/class.html

                    CFD高性能服务器 http://dyfluid.com/servers.html

                    深蓝 1 条回复 最后回复 回复 引用
                    • 深蓝
                      深蓝 @李东岳 最后由 编辑

                      @李东岳 是的,确实是这样的

                      bestucan 1 条回复 最后回复 回复 引用
                      • bestucan
                        bestucan 版主 副教授 @深蓝 最后由 编辑

                        @深蓝 你的问题已经超出我的能力范围了:zoule:
                        但是这里有个一个很完备的例子,应该对你有帮助。

                        一个简单而完整的UDS例子.zip

                        线性化,我能想到的方法就是泰勒级数展开了,不过要求值的变化范围不能太大,因为泰勒展开有个展开点,离展开点越近越准。这个过程应该是在软件外进行的,推出具体的计算方程再往UDF里写。

                        只能帮你到这了:xinxin2:

                        滚来滚去……~(~o ̄▽ ̄)~o 滚来滚去都不能让大家看出来我不是老师么 O_o

                        异步沟通方式(《posting style》from wiki)(下载后打开):
                        https://www.jianguoyun.com/p/Dc52X2sQsLv2BRiqnKYD
                        提问的智慧(github在gitee的镜像):
                        https://gitee.com/bestucan/How-To-Ask-Questions-The-Smart-Way/blob/master/README-zh_CN.md

                        深蓝 1 条回复 最后回复 回复 引用
                        • 深蓝
                          深蓝 @bestucan 最后由 编辑

                          @bestucan 谢谢您的回复
                          我已经找到解决办法了,在fluent中通过连续用C_UDSI(c,t,0)_G的办法不可行,又找到另外一种办法,如下:
                          未命名1631102135.png
                          是通过计算每个cell与周围cell的“浓度差”进行计算的。但是因为fluent只能获取到一个面左右两个网格(c0,c1)的信息,所以这种方法二阶导数最高只有二阶截断误差。想要有四阶阶段误差需要有9个节点的信息(二维),fluent的udf只能获得5个节点的信息。:ok2:

                          1 条回复 最后回复 回复 引用
                          • First post
                            Last post