CFD中文网

此问题同类于湍流动能耗散率的定义。湍流动能耗散率的定义为：
\begin{equation}
\varepsilon=2\nu_t\overline{S_{ij}S_{ij}}
\end{equation}
经过各向同性假定后有：
\begin{equation}
\varepsilon=\nu_t\overline{\frac{\p u_i'}{\p x_k}\frac{\p u_i'}{\p x_k}}
\end{equation}

@youmengtian 您好！我完整的格式是这样的:laile:

/*--------------------------------*- C++ -*----------------------------------*\ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org \\ / A nd | Version: 6 \\/ M anipulation | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; location "system"; object fvSchemes; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // fluxScheme Tadmor;// Kurganov; ddtSchemes { default Euler; } gradSchemes { default cellMDLimited Gauss linear 1.0; } divSchemes { default none; div(tauMC) Gauss linear; div(phi,nuTilda) Gauss linearUpwind limited; } laplacianSchemes { default Gauss linear corrected; } interpolationSchemes { default linear; } snGradSchemes { default corrected; } wallDist { method meshWave; }

感谢，好人一生平安

@东岳 这个边界条件已经整合到exteranlHeatFluxTemperature里面了，算的也是温度梯度，也就是$q=-\kappa \frac{\partial T}{\partial n}$，这种边界按照FDS技术手册里面应该是DNS的计算方法。
LES计算时给定的方法是$q=h \Delta T$,$h$是对流换热系数。openfoam里面的边界条件要么是固定值要么是梯度，这种差值不知道在of里面能不能实现

你好👋

有没有大佬做过电解碱水制氢的算例模拟啊，帮助文档里的电解水制氢复现了，点解碱水大佬们有没有做成功的case文件

@李东岳 李老师，就是感觉有时候，不能满足通过理论计算得到的边界层信息，比如不能满足边界层厚度啊，层数也达不到，这种一般怎么处理啊？感觉workbench的mesh虽然可以，但是我希望通过这个方式画网格，这样对我后续的计算方便一点。💋

@Dingcy 我再推荐几本书哈，一本Patankar的Numerical Heat Transfer and Fluid Flow，印象中前半本针对finite difference，后半本针对finite volume（这本特别适合初学理论又不想看特别难的，相当平实的语言）；还有一本Tannehill, Anderson, Pletcher的Computational Fluid Mechanics and Heat Transfer，这本可以当工具书用，大而全，倒不适合使劲的啃。当然还有Ferziger, Peric的那本Computational Methods for Fluid Dynamics，这本我个人感觉挺好的，就是需要你对张量符号、运算之类的比较熟悉，不然看起来经常会纠结于细节上的推导而忽略了物理和数学过程。

你用了两年半OpenFOAM了，应该是想让它实现什么功能都能很快做到（或者至少知道怎么去找相关的资料），不如在数值计算的理论方面多下下功夫。个人建议仅供参考。

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

@Foamer24

你回贴的这俩位我还有点印象，2016年的时候他们俩研究的还很深刻，后来就没怎么来过。不知道最近怎么样了，有没有在CFD。我有种预感他俩是不是转行了...

@东岳 万分感谢!

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

@bestucan 谢谢您！学习了

求助各位老师，在扩展出口计算域计算时，发现了如下问题，具体细节如下图，不同方案中红色截面处的总压大小不同，是什么原因造成的。（尤其是方案2和方案3总压的不同是因为什么原因）
问.png

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

微信截图_20201227163302.png
如图，这个代码是怎么实现的，原理是什么

https://www.bilibili.com/video/BV1or421x7Ut/?spm_id_from=333.999.0.0&vd_source=5dfca494f85875fa7f031a9f632b1a89
第9分钟介绍这种High Aspect Ratio对于瞬态、稳态的具体影响。

我觉得不能。RANS大部分都抹掉了。
有大佬看到系统的论述么？

@李东岳 谢谢东哥，就是我把resolveFeatrueAngle调大之后画网格没有捕捉这条特征线，画几何的时候因为在两个不同平面上，这条特征线是存在的

@liailei0627 在 Openfoam输出文件 中说：

有没有办法直接编程非常简单的，ASCII，没有括号，和上面的字头等，只是单纯的数字。

太简单了。你把openfoam当成C或者C++，在C里面怎么写代码就怎么写啊。

只要你不用openfoam语言输出（例如auto_write这种），完全可行。