CFD中文网

各位大佬，请问CFD-DEM coupling时，流体的控制方程组求解，用耦合式解法还是分离式解法？
非常感谢！

整体划分，然后删除多余的block，看看好用不

@Shane 非常感谢

这个啊，先排错，
看看是“只有20可以”，还是“数字小了才可以”。
比如试试19，21，试试200。

错误定的越准，描述才越准，也越容易搜到前人的经验

@Alex 太厉害了，大神，我想学习icem，您还有别的类似的icem实例教程吗，不胜感激

列表

@李东岳 这点我清楚，实际上这么多参数里也就总摩擦系数对不上，其他的基本上全对上了。LES的话得在文献里是说在入口用vortex method施加扰动，而且模拟的是有肋的，所以基本上不担心涡的问题。流向周期性的模型也可以在RANS模拟之后用TUI命令施加扰动，之前试过了涡能长期存在，不过那个case没算完结果对不对就不清楚了。

@HydroJames 我暂时还查不到原因，测试过程中我还引入了其他差异，比如我进行过换硬盘的操作

@Leon-Zhao :xiexie: 一键三连，哈哈哈哈，期待某天能在b站恰饭

采用overset重叠网格和UDF函数模拟串列三圆柱的涡激振动，其中UDF分别尝试了Newmark-Beta方法和4阶Runge Kutta法获取圆柱的振动响应，通过DEFINE_CG_MOTION宏赋予三个圆柱及component cells的运动速度，算例的雷诺数约200，采用k-omega sst模型，考虑水作为来流介质，计算过程中尝试了时间步长从1.0e-3缩小到1.0e-5等多个量级，但是求解过程总是出现升力和升力矩突然骤增，继而导致圆柱运动速度过大，最终计算发散。

请各位大佬帮忙看看是哪里出问题？

具体的UDF和部分设置如下：

#include "udf.h" #include "sg_mem.h" #include "dynamesh_tools.h" #define PI 3.141592654 #define zoneID_1 4 #define zoneID_2 16 #define zoneID_3 20 FILE *outNB,*outRK; static real y = 0.0; static real yRK = 0.0; static real dy = 0.0; static real vy = 0.0; static real vyRK = 0.0; static real vyRK2 = 0.0; static real ay = 0.0; static real current_time = 5; static real y2 = 0.0; static real y2RK = 0.0; static real dy2 = 0.0; static real vy2 = 0.0; static real vy2RK = 0.0; static real vy2RK2 = 0.0; static real ay2 = 0.0; static real current_time2 = 5; static real y3 = 0.0; static real y3RK = 0.0; static real dy3 = 0.0; static real vy3 = 0.0; static real vy3RK = 0.0; static real vy3RK2 = 0.0; static real ay3 = 0.0; static real current_time3 = 5; DEFINE_CG_MOTION(cylinder_1,dt,vel,omega,time,dtime) { real ctime = RP_Get_Real("flow-time"); real ctimestep = RP_Get_Integer("time-step"); real niter = N_ITER; if (current_time < ctimestep) { current_time = ctimestep; /*Define variables*/ /*Mesh variables*/ real cg[3],vcg[3]; /*Cylinder variables*/ real diameter = 0.063; real fn = 1.0892; real density = 998.2; real length = 1; real water_depth = 1; real mass_ratio = 0.3937; real damping_ratio = 0.01; real mass = mass_ratio*density*pow((0.5*diameter),2)*PI*length; real ad_mass = mass*(0); /*density*pow((0.5*diameter),2)*PI*water_depth;*/ real total_mass = mass + ad_mass; real k = 4*pow((PI*fn),2)*total_mass; real c = 2 * damping_ratio * sqrt(k*total_mass); /*Force calculation. Force = F_pressure + F_viscous*/ real fy = 0.0; real fvy = 0.0; int i; #if !RP_HOST Thread *tc,*thread; Domain *d = Get_Domain(1); face_t f; tc = Lookup_Thread(d,zoneID_1); thread = DT_THREAD(dt); NV_S(vel, =, 0.0); NV_S(omega, =, 0.0); real NV_VEC(A); begin_f_loop(f,tc) { if (PRINCIPAL_FACE_P(f,tc)) { fvy = F_STORAGE_R_N3V(f,tc,SV_WALL_SHEAR)[1]*-1; /*“*-1”表示方向*/ F_AREA(A,f,tc); /*Force calculation with a depth of 1m*/ fy += F_P(f,tc)*A[1] + fvy; } } end_f_loop(f,tc) #endif #if RP_NODE fy = PRF_GRSUM1(fy); #endif /*Dynamic mesh position*/ #if!RP_HOST for (i=0;i<3;i++) { cg[i]=DT_CG(dt)[i]; vcg[i] = DT_VEL_CG(dt)[i]; } Message("Position CG: %f \n",cg[1]); #endif node_to_host_real_2(fy,cg[1]); /*Numerical methods*/ /*Numark-beta*/ real beta = 0.25; real gamma = 0.5; real term0 = (1/(beta*dtime*dtime))*(mass+ad_mass) + (gamma/(beta*dtime))*c; real term1 = (1/(beta*dtime))*(mass+ad_mass) + ((gamma/beta)-1)*c; real term2 = ((1/(2*beta))-1)*(mass+ad_mass) + dtime*((gamma/(2*beta))-1)*c; real Keff = k + term0; real Reff = fy*water_depth + term0*cg[1] + term1*vy + term2*ay; Message("Velocity: %f \n",vy); dy = Reff/Keff - cg[1]; y += dy; real vprev = vy; vy = (gamma/(beta*dtime))*dy + (1-(gamma/beta))*vy + dtime*(1-(gamma/(2*beta)))*ay; ay = (1/(beta*dtime*dtime))*dy - (1/(beta*dtime))*vprev - ((1/(2*beta))-1)*ay; /*Runge-kutta 4th order*/ real K1 = (fy*water_depth - c*vyRK - k*yRK) / total_mass; real K2 = (fy*water_depth - c*(vyRK+dtime*0.5*K1) - k*(yRK+dtime*0.5*vyRK)) / total_mass; real K3 = (fy*water_depth - c*(vyRK+dtime*0.5*K2) - k*(yRK+dtime*0.5*vyRK+dtime*dtime*K1/4)) / total_mass; real K4 = (fy*water_depth - c*(vyRK+dtime*K3) - k*(yRK+dtime*vyRK+dtime*dtime*K1/2)) / total_mass; yRK = yRK + vyRK*dtime + dtime*dtime*(K1 + K2 + K3 + K4)/6; vyRK = vyRK + dtime*(K1 + K2 + K3 + K4)/6; /*Transfer result to the dynamic mesh*/ vel[0] = 0.0; vel[1] = vyRK; /*Save files*/ #if !RP_NODE /*Message ("Force = %f, pos = %f, vel = %f, acc = %f\n", fy, cg[1], y, vy);*/ if(NULL == (outNB = fopen("dataNB1.txt","a"))) { Error("Could not open file for append!\n"); } fprintf(outNB,"%16.4e %12.1f %16.3e %16.7f %16.7f %16.7f \n", ctime,niter, fy , cg[1], y, vy); fclose(outNB); if(NULL == (outRK = fopen("dataRK1.txt","a"))) { Error("Could not open file for append!\n"); } fprintf(outRK,"%16.4e %12.1f %16.3e %16.7f %16.7f %16.7f \n", ctime,niter, fy , cg[1], yRK, vyRK); fclose(outRK); #endif } /*Transfer result to the dynamic mesh*/ vel[0] = 0.0; vel[1] = vyRK; } DEFINE_CG_MOTION(cylinder_1_frontgrid_1,dt,vel,omega,time,dtime) { NV_S(vel, =, 0.0); NV_S(omega, =, 0.0); vel[0]=0.0; vel[1]=vyRK; } DEFINE_CG_MOTION(cylinder_1_overset_2,dt,vel,omega,time,dtime) { NV_S(vel, =, 0.0); NV_S(omega, =, 0.0); vel[0]=0.0; vel[1]=vyRK; } DEFINE_ZONE_MOTION(cylinder_1_zone,omega,axis,origin,velocity,time,dtime) { N3V_D(velocity, =, 0, 0, 0); N3V_S(origin, =, -0.32); N3V_D(axis, =, 0.0, 0.0, 1.0); velocity[1]=vyRK; } DEFINE_CG_MOTION(cylinder_2,dt,vel,omega,time,dtime) { real ctime = RP_Get_Real("flow-time"); real ctimestep = RP_Get_Integer("time-step"); real niter = N_ITER; if (current_time2 < ctimestep) { current_time2 = ctimestep; /*Define variables*/ /*Mesh variables*/ real cg[3],vcg[3]; /*Cylinder variables*/ real diameter = 0.063; real fn = 1.0892; real density = 998.2; real length = 1; real water_depth = 1; real mass_ratio = 0.3937; real damping_ratio = 0.01; real mass = mass_ratio*density*pow((0.5*diameter),2)*PI*length; real ad_mass = mass*(0); /*density*pow((0.5*diameter),2)*PI*water_depth;*/ real total_mass = mass + ad_mass; real k = 4*pow((PI*fn),2)*total_mass; real c = 2 * damping_ratio * sqrt(k*total_mass); /*Force calculation. Force = F_pressure + F_viscous*/ real fy = 0.0; real fvy = 0.0; int i; #if !RP_HOST Thread *tc,*thread; Domain *d = Get_Domain(1); face_t f; tc = Lookup_Thread(d,zoneID_2); thread = DT_THREAD(dt); NV_S(vel, =, 0.0); NV_S(omega, =, 0.0); real NV_VEC(A); begin_f_loop(f,tc) { if (PRINCIPAL_FACE_P(f,tc)) { fvy = F_STORAGE_R_N3V(f,tc,SV_WALL_SHEAR)[1]*-1; /*“*-1”表示方向*/ F_AREA(A,f,tc); /*Force calculation with a depth of 1m*/ fy += F_P(f,tc)*A[1] + fvy; } } end_f_loop(f,tc) #endif #if RP_NODE fy = PRF_GRSUM1(fy); #endif /*Dynamic mesh position*/ #if!RP_HOST for (i=0;i<3;i++) { cg[i]=DT_CG(dt)[i]; vcg[i] = DT_VEL_CG(dt)[i]; } Message("Position CG: %f \n",cg[1]); #endif node_to_host_real_2(fy,cg[1]); /*Numerical methods*/ /*Numark-beta*/ real beta = 0.25; real gamma = 0.5; real term0 = (1/(beta*dtime*dtime))*(mass+ad_mass) + (gamma/(beta*dtime))*c; real term1 = (1/(beta*dtime))*(mass+ad_mass) + ((gamma/beta)-1)*c; real term2 = ((1/(2*beta))-1)*(mass+ad_mass) + dtime*((gamma/(2*beta))-1)*c; real Keff = k + term0; real Reff = fy*water_depth + term0*cg[1] + term1*vy2 + term2*ay2; Message("Velocity: %f \n",vy2); dy2 = Reff/Keff - cg[1]; y2 += dy2; real vprev = vy2; vy2 = (gamma/(beta*dtime))*dy2 + (1-(gamma/beta))*vy2 + dtime*(1-(gamma/(2*beta)))*ay2; ay2 = (1/(beta*dtime*dtime))*dy2 - (1/(beta*dtime))*vprev - ((1/(2*beta))-1)*ay2; /*Runge-kutta 4th order*/ real K1 = (fy*water_depth - c*vy2RK - k*yRK) / total_mass; real K2 = (fy*water_depth - c*(vy2RK+dtime*0.5*K1) - k*(y2RK+dtime*0.5*vy2RK)) / total_mass; real K3 = (fy*water_depth - c*(vy2RK+dtime*0.5*K2) - k*(y2RK+dtime*0.5*vy2RK+dtime*dtime*K1/4)) / total_mass; real K4 = (fy*water_depth - c*(vy2RK+dtime*K3) - k*(y2RK+dtime*vy2RK+dtime*dtime*K1/2)) / total_mass; y2RK = y2RK + vy2RK*dtime + dtime*dtime*(K1 + K2 + K3 + K4)/6; vy2RK = vy2RK + dtime*(K1 + K2 + K3 + K4)/6; /*Transfer result to the dynamic mesh*/ vel[0] = 0.0; vel[1] = vy2RK; /*Save files*/ #if !RP_NODE /*Message ("Force = %f, pos = %f, vel = %f, acc = %f\n", fy, cg[1], y, vy);*/ if(NULL == (outNB = fopen("dataNB2.txt","a"))) { Error("Could not open file for append!\n"); } fprintf(outNB,"%16.4e %12.1f %16.3e %16.7f %16.7f %16.7f \n", ctime,niter, fy , cg[1], y2, vy2); fclose(outNB); if(NULL == (outRK = fopen("dataRK2.txt","a"))) { Error("Could not open file for append!\n"); } fprintf(outRK,"%16.4e %12.1f %16.3e %16.7f %16.7f %16.7f \n", ctime,niter, fy , cg[1], y2RK, vy2RK); fclose(outRK); #endif } /*Transfer result to the dynamic mesh*/ vel[0] = 0.0; vel[1] = vy2RK; } DEFINE_CG_MOTION(cylinder_2_frontgrid_1,dt,vel,omega,time,dtime) { NV_S(vel, =, 0.0); NV_S(omega, =, 0.0); vel[0]=0.0; vel[1]=vy2RK; } DEFINE_CG_MOTION(cylinder_2_overset_2,dt,vel,omega,time,dtime) { NV_S(vel, =, 0.0); NV_S(omega, =, 0.0); vel[0]=0.0; vel[1]=vy2RK; } DEFINE_ZONE_MOTION(cylinder_2_zone,omega,axis,origin,velocity,time,dtime) { N3V_D(velocity, =, 0, 0, 0); N3V_S(origin, =, 0.0); N3V_D(axis, =, 0.0, 0.0, 1.0); velocity[1]=vy2RK; } DEFINE_CG_MOTION(cylinder_3,dt,vel,omega,time,dtime) { real ctime = RP_Get_Real("flow-time"); real ctimestep = RP_Get_Integer("time-step"); real niter = N_ITER; if (current_time3 < ctimestep) { current_time3 = ctimestep; /*Define variables*/ /*Mesh variables*/ real cg[3],vcg[3]; /*Cylinder variables*/ real diameter = 0.063; real fn = 1.0892; real density = 998.2; real length = 1; real water_depth = 1; real mass_ratio = 0.3937; real damping_ratio = 0.01; real mass = mass_ratio*density*pow((0.5*diameter),2)*PI*length; real ad_mass = mass*(0); /*density*pow((0.5*diameter),2)*PI*water_depth;*/ real total_mass = mass + ad_mass; real k = 4*pow((PI*fn),2)*total_mass; real c = 2 * damping_ratio * sqrt(k*total_mass); /*Force calculation. Force = F_pressure + F_viscous*/ real fy = 0.0; real fvy = 0.0; int i; #if !RP_HOST Thread *tc,*thread; Domain *d = Get_Domain(1); face_t f; tc = Lookup_Thread(d,zoneID_3); thread = DT_THREAD(dt); NV_S(vel, =, 0.0); NV_S(omega, =, 0.0); real NV_VEC(A); begin_f_loop(f,tc) { if (PRINCIPAL_FACE_P(f,tc)) { fvy = F_STORAGE_R_N3V(f,tc,SV_WALL_SHEAR)[1]*-1; /*“*-1”表示方向*/ F_AREA(A,f,tc); /*Force calculation with a depth of 1m*/ fy += F_P(f,tc)*A[1] + fvy; } } end_f_loop(f,tc) #endif #if RP_NODE fy = PRF_GRSUM1(fy); #endif /*Dynamic mesh position*/ #if!RP_HOST for (i=0;i<3;i++) { cg[i]=DT_CG(dt)[i]; vcg[i] = DT_VEL_CG(dt)[i]; } Message("Position CG: %f \n",cg[1]); #endif node_to_host_real_2(fy,cg[1]); /*Numerical methods*/ /*Numark-beta*/ real beta = 0.25; real gamma = 0.5; real term0 = (1/(beta*dtime*dtime))*(mass+ad_mass) + (gamma/(beta*dtime))*c; real term1 = (1/(beta*dtime))*(mass+ad_mass) + ((gamma/beta)-1)*c; real term2 = ((1/(2*beta))-1)*(mass+ad_mass) + dtime*((gamma/(2*beta))-1)*c; real Keff = k + term0; real Reff = fy*water_depth + term0*cg[1] + term1*vy3 + term2*ay3; Message("Velocity: %f \n",vy3); dy3 = Reff/Keff - cg[1]; y3 += dy3; real vprev = vy3; vy3 = (gamma/(beta*dtime))*dy3 + (1-(gamma/beta))*vy3 + dtime*(1-(gamma/(2*beta)))*ay3; ay3 = (1/(beta*dtime*dtime))*dy3 - (1/(beta*dtime))*vprev - ((1/(2*beta))-1)*ay3; /*Runge-kutta 4th order*/ real K1 = (fy*water_depth - c*vy3RK - k*y3RK) / total_mass; real K2 = (fy*water_depth - c*(vy3RK+dtime*0.5*K1) - k*(y3RK+dtime*0.5*vy3RK)) / total_mass; real K3 = (fy*water_depth - c*(vy3RK+dtime*0.5*K2) - k*(y3RK+dtime*0.5*vy3RK+dtime*dtime*K1/4)) / total_mass; real K4 = (fy*water_depth - c*(vy3RK+dtime*K3) - k*(y3RK+dtime*vy3RK+dtime*dtime*K1/2)) / total_mass; y3RK = y3RK + vy3RK*dtime + dtime*dtime*(K1 + K2 + K3 + K4)/6; vy3RK = vy3RK + dtime*(K1 + K2 + K3 + K4)/6; /*Transfer result to the dynamic mesh*/ vel[0] = 0.0; vel[1] = vy3RK; /*Save files*/ #if !RP_NODE /*Message ("Force = %f, pos = %f, vel = %f, acc = %f\n", fy, cg[1], y, vy);*/ if(NULL == (outNB = fopen("dataNB3.txt","a"))) { Error("Could not open file for append!\n"); } fprintf(outNB,"%16.4e %12.1f %16.3e %16.7f %16.7f %16.7f \n", ctime,niter, fy , cg[1], y3, vy3); fclose(outNB); if(NULL == (outRK = fopen("dataRK3.txt","a"))) { Error("Could not open file for append!\n"); } fprintf(outRK,"%16.4e %12.1f %16.3e %16.7f %16.7f %16.7f \n", ctime,niter, fy , cg[1], y3RK, vy3RK); fclose(outRK); #endif } /*Transfer result to the dynamic mesh*/ vel[0] = 0.0; vel[1] = vy3RK; } DEFINE_CG_MOTION(cylinder_3_frontgrid_1,dt,vel,omega,time,dtime) { NV_S(vel, =, 0.0); NV_S(omega, =, 0.0); vel[0]=0.0; vel[1]=vy3RK; } DEFINE_CG_MOTION(cylinder_3_overset_2,dt,vel,omega,time,dtime) { NV_S(vel, =, 0.0); NV_S(omega, =, 0.0); vel[0]=0.0; vel[1]=vy3RK; } DEFINE_ZONE_MOTION(cylinder_3_zone,omega,axis,origin,velocity,time,dtime) { N3V_D(velocity, =, 0, 0, 0); N3V_S(origin, =, 0.32); N3V_D(axis, =, 0.0, 0.0, 1.0); velocity[1]=vy3RK; }

运动速度.png
运动速度

运动位移.png
运动位移

压力系数.png
压力系数

动网格设置.png
动网格设置

#问题1 自我解决，感觉像是bug.
0_1488141087886_streamline.png

个人觉得这个功能特别实用，特别是对于边界很多的案例，能直观看到边界在哪个位置。但是2020以后的版本都没有找到此功能，想问一下是阉割掉了还是被别的功能替代了？谢谢！微信图片_20230818185948.png

简单的点线面操作撤销是没问题的，检查拓扑或者是操作包含的数量太多ICEM就没办法了，解决这个问题的最好办法就是你在做包含多步操作的问题时候先另存为一下是最保险的，问了一下学长他也是用这个方法的

好像是分内部点（所有面都不是边界）和边界点，内部点根据cell插值（距离倒数分之一），边界点根据face插值。

@晴空星茗 我不是老师，我也是研究僧:xinxin3: 。“做正确的事就有好的回报”的环境不是白来的。你想要这样的环境，别人也想要这样的环境，可是得有人去创造。白白享受这样的环境无异于蚕食巨人的血肉。努力创造这样的环境，无论在多大、多小的范围内。己欲立而立人，己欲达而达人。不信你看本论坛的“青年成长计划”、“验证算例支持”:shangxue:

对待那些为难你的形势最好的报复就是，当我们走到了那个位置还记得现在的心情和要做什么。

各位CFD大佬，小弟最近在用fluent做模拟时遇到一个问题，物理背景就是管内的热水以射流的方式喷出然后被周围的环境冷水稀释，在做小模型模拟（x=40m y=2m z=13.3cm）的时候，模拟结果尚可，对称性很完美，但在做大模型的时候（x=4000m y=2000m z=10m）,模拟结果很不好。
用的稳态模拟，大模型网格数量大概七百多万，用的reliazable k-e湍流模型，进口的湍流参数设置用CFD计算器估算过了，另外在模拟的时候fluent一直有湍流粘度比超限的提示，但小模型就不会有，不知道最后的不对称结果和这个是不是有关系。
1564dc71-8ed2-4905-ba1b-67fbca1936e3-图片.png
热水从圆的四周均匀喷射出来，然后被流动的冷水稀释，冷水流动方向从左至右，可以看到图中并没有达到一个对称效果，等值线都没有对称。
21f863e0-2b3c-4e9f-994a-ae3797889acf-图片.png
热水进口在上图的最底部，垂直向上喷射，中间比较黑的那一块有一块顶盖，作用是让热水向周围扩散开来而不是垂直向上冲，图中看到网格有些瑕疵的地方是tecplot的显示问题，可以忽略。
a652dbd8-6ce0-4c30-9c4a-97e6fee8e965-892fea64760b5fcdfe9299ba371783e.png
从这张图可以看出热水向上喷射，被顶盖挡住，最后向四周扩散上浮的过程。

总结一下几点疑问：
1.模拟时fluent的湍流粘度比超限警告是否和最后模拟结果的不对称有关系？
2.小模型模拟结果对称，而大模型就不行，和网格的大小有关系吗？小模型的最大网格大概也就在厘米级别，而大模型的最大网格一般在米级别（最大有五米左右）。
其实最核心的诉求还是想解决对称性的问题，最好不要用对称边界条件，能用整个模型就用整个模型。

是的，理想气体，明天我更新下。

另外，温度本身不是守恒变量，因此，CFD一般求解守恒变量，但是温度方程通常比较简单，一些不太成熟的代码会使用温度防尘

:chouchou:

之前也困惑问过类似的问题，李老师说复杂的情况用ICEM 之后网格转换一下，这样确实是省时又省力哦。

@cfd-china 谢谢回复，我研究了源码已经搞懂了，顺便附上我总结的链接laplacian项离散解码