CFD中文网

@队长别开枪 感谢感谢 😊

如题，这段时间在做一个分子模拟的案例，需要应用到N-S方程反求粘度，想知道这里边的一些基本思路。

0_1515741447788_草图.png

小的时候我妈也跟我说过。我家条件不好，但是他们给我创造了一个安心的学习环境。那个时候很不理解这句话：

“家里什么事都不需要你操心，抓紧时间学习吧”

直到30而立，才明白，有时间学习，或者说有时间、有钱去支持你从事科研，非常难得。

小时候我们是在学习，支撑学习的经费是家里给的。现在我们也是在学习，科研何尝不是学习呢？但是我们需要在学习的同时，自己去寻找支撑自己学习的经费。

现在的高中生或者本科生们，珍惜时间好好学习吧。

@bestucan 在 关于OpenFOAM计算的是否是动压力系数？ 中说：

都不可压了，计算中密度也不会变。可以看看东岳流体对基本求解器的解析。

密度是设定的，在 constant/transportProperties 里

谢谢回复！

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

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

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

As an editor of academic papers and articles, I would interpret “intensively” in your ultimate citation as indicating that these nanostructures have been studied in depth, whereas “extensively” informs us that numerous studies were undertaken during the decade in question. In other words, “intensively” reflects degree and “extensively” frequency. Thanks for an informative exposition, I look forward to reading more.

二变量高斯分布：
\begin{equation}
f(u,v)=\frac{1}{2\pi\sigma_1\sigma_2\sqrt{1-\rho^2}}\exp\left(-\frac{1}{2(1-\rho^2)}\left[\frac{(u-\mu_1)^2}{\sigma_1^2}-\frac{2\rho(u-\mu_1)(v-\mu_2)}{\sigma_1\sigma_2}+\frac{(v-\mu_2)^2}{\sigma_2^2}\right]\right)
\end{equation}
MGF为：
\begin{equation}
m_{i,j}=\exp\left(i\mu_1+j\mu_2+0.5(\sigma_1^2i^2+2\rho\sigma_1\sigma_2ij+\sigma_2^2j^2)\right)
\end{equation}
纯矩计算方法为
\begin{equation}
\begin{split}
m_{0,0}&=1\\
m_{1,0}&=\mu\\
m_{2,0}&=\mu^2+\sigma^2\\
m_{3,0}&=\mu^3+3\mu\sigma^2\\
\end{split}
\end{equation}
假设$\mu_1=10,\mu_2=20,\sigma_1=\sigma_2=2,\rho=0.5$，有纯矩：
\begin{split}
m_{0,0}&=1\\
m_{1,0}&=10\\
m_{2,0}&=104\\
m_{3,0}&=1120\\
m_{0,1}&=20\\
m_{0,2}&=404\\
m_{0,3}&=8240\\
\end{split}
同时有混合矩
\begin{equation}
m_{i,j}=\exp\left(36\right)
\end{equation}
不行，混合矩计算方法不对

meanVelocityForce只用于不可压缩，但是涉及传热又不区分可压缩与不可压缩，那么像buoyantPimpleFoam这种求解器可以用meanVelocityForce来驱动么？有了解的前辈解释下么？

patchInternalField表示你说的这些量

@bestucan 嗯嗯，感谢大佬。检查过，check和Quallity其实都没有问题，再次感谢~

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

考虑最简单的动量交换，液相：
\begin{equation}
\frac{\rd \bfU_\rc}{\rd t}=-\bfA
\end{equation}
气相：
\begin{equation}
\frac{\rd \bfU_\rd}{\rd t}=\bfA
\end{equation}
乍一看，如果$\bfA$是负的，那么会导致$\bfU_\rd$向下走（为负），$\bfU_\rc$（为正）。但是这并不符合物理，考虑一个管子的颗粒，如果颗粒向下走，必然会带动空气同时向下走。但是方程缺不是这种体现，为什么呢？

实际上，上述方程并没有写完整，完整形式应该是这样：
液相：
\begin{equation}
\frac{\rd \bfU_\rc}{\rd t}=-\bfA=-\frac{1}{\tau}(\bfU_\rc-\bfU_\rd)
\end{equation}
气相：
\begin{equation}
\frac{\rd \bfU_\rd}{\rd t}=\bfA=\frac{1}{\tau}(\bfU_\rc-\bfU_\rd)
\end{equation}
如果开始的时候$\bfU_\rc=0$，则为

12岁，考个好初中

15岁，考个好高中

18岁，考个好大学

22岁，考研

25岁，找个好工作。

27岁，相亲组建个家庭

30岁，挣钱攒钱买房 -- 我目前在这里 :surprised: 10万一平 :upset:

40岁，养孩子养老人，养车养房。

50岁，给孩子买房子

60岁，带孙子养孙子

70岁，是不是可以享福了

@李东岳 好的，谢谢李老师耐心解答

你需要系统看一下流体力学的相关知识，可能你没有这方面的背景所以不太理解。对物理空间中的一块区域而言，其所包含的物理量有几种变化方式，要么它自己随时间在变化（变化率，比如密度因为温度升高降低了），要么因为有流动带它进来或者出去（对流，比如能量），要么它自己在向外扩散或者外部在向这块区域扩散（扩散，比如这一块是高浓度区域而周围是低浓度），要么这块区域中这个物理量自己在生成或者消灭（源项，比如有化学反应）。所有的这些变化，都需要满足物理学的基础即守恒定律，其数学表达就是所谓的“广义”NS方程。这里“广义”是指这个方程描述了物理量的一般过程。

phi等于1，这个方程代表质量守恒；phi等于速度，这个方程代表动量守恒；phi等于内能，这个方程代表能量守恒；phi等于某个物质的浓度，那么这个方程就代表该物质的质量守恒，等等。

链接文本

参考页面左边栏。

@东岳 谢谢李老师的回复！
我看了下跟公式（7）到（8）的过程很相似，但是为什么公式7到公式8的过程中将公式7左边第四项给略去呢？

@dawei 手动指定库文件路径后编译安装成功