CFD中文网

reFineMesh可以加密一个方向，这个需要指定区域，你试试topoSet+refineMesh呢》

https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1compressible_1_1humidityTemperatureCoupledMixedFvPatchScalarField.html

https://www.cfd-online.com/Forums/openfoam-solving/147129-condensation-using-openfoam.html

https://www.researchgate.net/publication/279643607_Modeling_of_Direct_Contact_Condensation_With_OpenFOAM

https://jhmtr.semnan.ac.ir/article_467.html

https://github.com/NimaSam/phaseChangeHeatFoam

wps的国际版更干净。但是国内版的缺失字体会自动下载，国际版好像因为版权什么的，这类功能少点。不知道老师用的哪版

@Yan 您好，请问您是用的什么后处理软件，在导出速度的时候可以把笛卡尔坐标系改成柱坐标系

“对流”-“扩散”型、偏导数型（对时间t的偏导数加对xyz的偏导数）

没太理解这个描述。

此外，湍流输运方程（如k方程、omega方程）中的dissipation term中也不包含偏导数

在源项里面

采用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
动网格设置

说明前期没收敛啊，网格无关性是要做的

@卡洛 请问你是用的哪种方法？

https://www.cfd-online.com/Forums/openfoam-meshing/133850-snappyhexmesh-problems.html

截断误差首项阶次一样

我推了一下，CD的更小。不一样啊。

捕获.PNG

你是不是想把曲面展开为一个平面？https://blog.csdn.net/mrbaolong/article/details/106949340

@Junren-Hou 在 ICEM划分网格 中说：

可以试试check block，有些时候其实没有问题，就它软件自己的毛病，然后block里最后有一个check block，可以让他自己检查一下。里面有好多选项，都试试，什么check/fixed block

回忆了一下，我的旋转区域用ICEM画的应该没有问题。出问题的应该是我画的外流场区域，一方面追求地面边界层，致使网格尺度过小，一方面为了平衡网格数量，过渡不均匀。猜测应该是外流区域的地面边界层设置导致了相关报错的产生。

可能是网格太大，内存不够用。换个文件尺寸小点的网格文件试试

topo.png
我需要用ICEM画multiblock网格，几何类似于上图的模型，我知道三角形Y-block和圆形O-block，那请问这两个如何结合起来呢？在这种结构中它的edge并不共用相同的vertices
第二个问题是icem的multiblock的网格如何光顺？我尝试使用了pre-mesh smooth，但是这个功能很慢，容易卡死。我只知道可以用edit edge调整spline来改变edge形状，如果一个一个edge调是不是有点笨？

@李东岳 十分感谢！

@李东岳 忘了说了，网格都是结构化的，长方体网格:duang:

我采用的稳态MRF模型进行模拟，如果桨距角改变，同一来流速度下给旋转域设置的旋转速度需要改变吗

各位大佬好，模型是一个0.125x0.5x0.6m的光滑通道，上下壁面均为无滑移的300K恒温壁面，进口温度320K，进出口及左右壁面均采用周期性边界条件，模拟步骤是先用RANS模拟，之后使用TUI指令添加扰动然后切换到LES模拟。按理说流体温度应该小于等于320K而且随时间逐渐减小，但是云图画出来最高温度有323K，又模拟了0.5s后最高温度没什么变化其他地方的温度似乎也没什么变化。周期性设置如图（压力梯度最开始为0）搜狗截图20230903172308.png
我想知道是为啥 :136:

很多湍流模型的假定就是各向同性，这样模化起来、封闭起来比较简单。各向异性那种，就得上雷诺应力模型了，明显要复杂

对称几何对称网格对称边界条件，但最后模拟出来的结果不对称的问题。

这个问题是会有发生，就类似圆柱扰流，注定是不对称的。不过如果是一个可以对称的物理，加上一个完全对称的数值设置。发生不对称的结果。就不太好说了。