CFD中文网

@羽之下 谢谢

@carmelosun 如果流体是理想气体之类，热物性是根据温度和压力计算得到，入口温度低，初始场温度高。其实入口速度还是0.4，但是由于入射的流体密度比初始场内密度大，就会有个加速效果，看起来像是入口速度大

@bestucan
感谢您的回复！“不同方向的流动计算差异就大”这个解释很直接明了，我突然就有点明白这个网格线出现的原因了:xinxin:

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

81a56a41-f2c6-4e5c-8504-2cb7df0a1f28-image.png ，我看帮助文档中DPM对droplet的描述是换热、蒸发和沸腾，这都是液滴材料从液相传质到气相，请问能不能实现从气相到液相的传质（就是吸收）？

@yuan_neu 建议看看《张量分析》的相关课程
二阶 【各向同性】的 【张量函数】 （ Dij(r) ）
只可能是上述 (6.25) 的表达形式

所以，你首先要明白“张量”，其次再搞清楚“张量函数”的定义，然后需要理解“各向同性”的概念，最后才能明白 6.25 式是怎么推导出来的

实际上，NS 方程的 Newton 本构关系，也是 “可以” 基于 【线性】和【各向同性】的假设直接写出来的 ！

推荐看看清华黄克智的《张量分析》教材

缺这个包：binutils

或者路径没配置对，可以用echo $PATH看看有没有/usr/bin这一项

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

@gsky0809 老哥，物质传输方程在哪里设置啊，我刚开始学宏观偏析，还不太懂，能指导一下吗

@冠竹 有可能是入口速度太小了

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

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

@tens 好的，谢谢你，我试一下

@dzw05 我是这么使用这个通量的，surfaceScalarField phiVp = fvc::flux(Vp); 计算时是fvm::div(phiVp, C)，然后使用求解器时要求我设置div(interpolate(Vp),C)的格式。如果使用surfaceScalarField phiVp = fvc::interpolate(Vp) & mesh.Sf();要求我设置div((interpolate(Vp)&S),C)的格式，所以我想问他们的区别。
谢谢老师

parcel。parcel守恒，然后改变particle数

@李东岳 coupled solver比segregated solver能更快收敛，但在处理复杂几何形状物体时不太容易收敛[1]。

[1] http://www.tfd.chalmers.se/~hani/kurser/OS_CFD_2012/KlasJareteg/KlasJareteg-OF2012-Project-3-Version-2.0.pdf

@东岳 感谢李老师！

weblate ，或者 github。
debian中文手册翻译
但学习成本都挺高。

但是特性好，版本控制啥的，weblate 就是个轻量的版本控制系统。po 文件可以统一术语。后台就是这个软件 gettex

曾经想参与 debian 手册的翻译，没时间入门这些东西。:xinlei:

也有会用软件一个人带着其他人翻，用邮件列表报翻好的片段，然后统一上传。

对这种不是那么重的翻译，用专业的方法成本太大。那种翻译框架搭好了是照着往多种语言翻译的。用土法对个人压力过大，这种压力不只是时间紧迫性和工作量上的。没有好的工具，hold 不住相当碎片化的翻译数据，邮件勉强可以，版本控制器最适合（全历史记录）

可以扣搜扣搜有没有日文版和日文翻译者，取取经:chitang: 一般好多技术文档有国际翻译，八成就有日语翻译。
日本 OpenFOAM 大本营：googlegroups
组织上决定派你去踩点，加油:chigua2:

网上实在搜不到
作者: J.O. Hinze (Author)
出版社: McGraw-Hill Inc.,US; 2nd Revised edition edition

ISBN: 9780070290372 出版时间: 1975-05

邮箱 yuan_neu@163.com

好像我问的有问题，$A=0.5(\overline{p_1}-\overline{p_0})\times(\overline{p_2}-\overline{p_0})$这个东西看起来是个矢量。但是面积不是矢量啊。我回去对对OpenFOAM的代码。

我看了下，OpenFOAM里面mesh.Sf()对应的是这个：$0.5(\overline{p_1}-\overline{p_0})\times(\overline{p_2}-\overline{p_0})$，mesh.magSf()对应的是这个：$\left|0.5(\overline{p_1}-\overline{p_0})\times(\overline{p_2}-\overline{p_0})\right|$，后面这个应该就是面积了。

@evensun $\left|0.5(\overline{p_1}-\overline{p_0})\times(\overline{p_2}-\overline{p_0})\right|$ 这个面积可以理解不