openfoam 不可压缩计算中的能量守恒
-
各位老师好，
我最近在用openfoam计算一个简单的强迫振动问题，想要验证能量守恒的关系(Loose Leaf for Modern Compressible Flow With Historical Perspective (John Anderson) )。
因为我是计算不可压缩的问题，所以能量方程有一定简化。这里我等号左侧第一二项为0（没有热量变化和旋转），左侧最后一项也为0（没有外力）
模型如下:
后面通过输出每个时刻的速度压力场处理公式里的各项，但是结果相差很大，可以说完全不守恒。我以为不可压缩问题天然满足能量守恒公式，结果和算例贴在下面。
c2case.7z这个结果相差一个数量级
我分析有两个原因，一个是选的能量公式 有误，在openfoam里输出每项的时候有误，所以在comsol里验证了一次，结果和openfoam的相差不大
（后面整理了再贴出来）。如何计算质量守恒和能量守恒 | COMSOL 博客感觉很奇怪，想请教各位老师。
-
@李东岳
李老师，这个项我理解是应力张量拆开了。
我计算的这项就是各个边界上的压力和速度通量。
另外不可压缩不存在这项的话是否意味着计算结果为0，但我算出来是有值的。
另外我参考comsol里的能量公式，去掉传热部分应该是一致的。
https://cn.comsol.com/blogs/how-to-calculate-mass-conservation-and-energy-balance -
我是输出各个边界上的速度压力场自己积分的这一项。代码里是B3,按照这样计算是有非零值的，不知道哪里出了问题。```
// 计算速度场的梯度 scalar rho = 1000; // 流体密度 const scalar mu = 1e-3; // 粘性系数 tmp<volTensorField> gradUTmp = fvc::grad(U); const volTensorField& gradU = gradUTmp(); // 初始化总的B3和剪切通量 scalar totalB3 = 0.0; scalar totalViscousFlux = 0.0; // 遍历所有边界patch forAll(mesh.boundaryMesh(), patchI) { const polyPatch& patch = mesh.boundaryMesh()[patchI]; word patchName = patch.name(); // 如果对特定的patch感兴趣，可以添加条件检查 if ((patchName == "front")||(patchName == "back")) { continue; // 跳过名为"front"或"back"的边界 } // 初始化当前patch的B3和剪切通量 scalar B3Patch = 0.0; scalar ViscousFluxPatch = 0.0; // 循环遍历patch上的所有面 forAll(patch, faceId) { const vector& Sf = patch.faceNormals()[faceId]; // 从patch对象获取归一化法向向量 scalar faceArea = mag(patch.faceAreas()[faceId]); // 计算面积向量的模长作为实际面积 const vector& U_boundary = U.boundaryField()[patchI][faceId]; const Tensor<double>& faceGradU = gradU.boundaryField()[patchI][faceId]; // 计算应力张量（包含法向和切向应力） Tensor<double> stressTensor = mu * (faceGradU + faceGradU.T()); // 计算壁面剪切应力（只保留切向分量） vector wallShearStress = stressTensor & Sf; wallShearStress -= (wallShearStress & Sf) * Sf / (magSqr(Sf) + VSMALL); // 壁面剪切应力与速度的点积 scalar wallShearStressDotU = wallShearStress & U_boundary; const volScalarField& pressure = mesh.lookupObject<volScalarField>("p"); scalar p_boundary = pressure.boundaryField()[patchI][faceId]; scalar pressureFlux = p_boundary * (U_boundary & Sf); // 根据面积加权的壁面剪切应力与速度的点积以及压力与速度的通量 ViscousFluxPatch += faceArea * wallShearStressDotU; B3Patch += faceArea * pressureFlux; } // 输出当前边界上的B3和剪切通量 Info << "Patch: " << patchName << ", Integrated B3: " << B3Patch << ", Integrated Viscous Flux: " << ViscousFluxPatch << endl; // 将当前patch的B3和剪切通量累加到总和中 totalB3 += B3Patch; totalViscousFlux += ViscousFluxPatch; } // 输出总的B3和剪切通量 Info << "Total Integrated B3 on all patches: " << totalB3 << endl; Info << "Total Integrated Viscous Flux on all patches: " << totalViscousFlux << endl; -
@李东岳 在 openfoam 不可压缩计算中的能量守恒 中说：
这个方程，如果看积分项，就是左边第4个，为什么存在压力的速度对流？不可压缩不存在这个
老师好，这个公式我又重新检查了，对应的是压力梯度项。
这部分展开后，在不可压缩流动中，右侧第一项确实是0，但是第二项积分应该不为零。同时我也很难理解，这两项积分后通过高斯公式转换为面积分不应该是相同的式子吗，怎么会有不一样的结果呢。
7/10