@李东岳
我是输出各个边界上的速度压力场自己积分的这一项。代码里是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;
