CFD中文网

1.根据现有的结果计算一下局部库朗数，看看是否是突变的速度导致的网格变形。（尝试调整一下时间项的离散方案）
2.用比较稀疏的网格尝试计算，看看能否承受变形（去掉边界层）
3.松弛因子accelerationRelaxation 0.4;可以再调整，网格变形方式也可以尝试调整。
大概是这些方向，如果都不行。那么可以考虑使用重叠网格方案，对于刚体运动是一定可行的。

@芋圆不圆

我对你这个方程不太熟悉，你自己可以对照原方程检查量纲。如果没弄错的话，
在 TEqn 方程中，fvm::ddt(T) 的维度应该是 [1 0 -3 1 0 0 0]（即温度的时间导数，单位是 K/s）。

fvm::div(phiFluid[i], T) 的维度应该是 [1 0 -3 1 0 0 0]（即对流项，单位是 K/s）。

fvm::laplacian(alphaEff, T) 的维度应该是 [1 0 -3 1 0 0 0]（即扩散项，单位是 K/s）。

rad.ST(rhoCp, T) 的维度应该是 [1 0 -3 1 0 0 0]（即辐射源项，单位是 K/s）

要让 rad.ST(rhoCp, T) 的返回值与方程中其他项的维度一致，可以检查 rad.ST 的实现或传入的参数(rhoCp)。大概是这个思路。

@xumengxin 对的，是这段时间窗口的平均值，对应的窗口值应该也有保存。

可以把方程各项的维度输出一下，有可能通量换算出问题了。
Info << "phiFluid[i] dimensions: " << phiFluid[i].dimensions() << endl;
Info << "alphaEff dimensions: " << alphaEff.dimensions() << endl;
Info << "rad.ST dimensions: " << rad.ST(rhoCp, T).dimensions() << endl;

用这个试一下

fieldAverage1
{
    type fieldAverage;
    functionObjectLibs ("libfieldFunctionObjects.so");
    resetOnRestart true;
    resetOnOutput false;
    outputControl   outputTime;
    fields
    (
        U
        {
            mean on;
            prime2Mean on;
            base time;//iteration
            windowType exact;//计算精确窗口平均值
            window 4000;
            allowRestart true;
        }
        p
        {
            mean on;
            prime2Mean on;
            base time;//iteration
            windowType exact;//计算精确窗口平均值
            window 4000;
            allowRestart true;
        }
    );
}

@Shihang-Chen
你好，我在做类似的代码。想请教伪时间步植入后，计算结果中出现松弛要怎么处理（计算的周期和实际周期不一致）。我是在动量方程中加入的隐式伪时间步。

@lwjetmann 多谢，我自己增加了一个solver，确实是可行的。

各位老师好，

我想要在运动系统中增加一个冲击位移，尝试了很多办法不能实现。
静止系统中的冲击可以用自定义边界条件来实现，我尝试了一个简单脉冲函数，可以很好的反映这个冲击。

如下在pointDisplacement文件中自定义待冲击的边界即可

{
        type            codedFixedValue;
        value           uniform (0 0 0); // 初始位移
        name            pulseDisplacement; // 自定义边界条件的名称
        code
        #{
        const scalar t = this->db().time().value();
        const scalar t0 = 5; // 位移脉冲的发生时刻
        const vector pulseAmplitude = vector(0, 0.001, 0); // 脉冲幅度

        if (t >= t0 && t < t0 + this->db().time().deltaT().value())
        {
            operator==(pulseAmplitude);
        }
        else
        {
            operator==(vector::zero);
        }
        #};
    }

但是6DOF运动中，pointDisplacement会被占用，他必须定义为sixdofsover里的边界（例如caculate）。这样上面静止系统的方案就不可行了。
我又尝试在pimpleFoam主程序中将位移写为附加项，想在6dof更新状态时加上去

if (currentTime >= pulseTime && currentTime < pulseTime + runTime.deltaT().value())
        {
            // 获取当前位移状态
            vector currentDisplacement = motion.centreOfRotation(); // 获取当前的中心位置

            // 将脉冲位移量叠加到当前位移
            vector newDisplacement = currentDisplacement + pulseDisplacement;

            // 使用 transform() 更新刚体的位置
            tensor identity(1, 0, 0, 0, 1, 0, 0, 0, 1);  // 单位张量，不涉及旋转
            // 调用合适的方法来更新位移
            motion.update(false, newDisplacement, vector(0, 0, 0), runTime.deltaT().value(), runTime.deltaT().value()); // 调用更新方法

			
            
			Info << "Pulse applied at time " << currentTime << " with displacement: " << pulseDisplacement << endl;
        }

这个方案也不可行，我估计是6dof不仅涉及到刚体位置的变换还有对应网格的更新。请各位老师帮忙看下。

各位老师好，

我想实现弹簧固定的模型上再耦合给定的运动，例如固定刚度的刚体在特定的时间范围上指定一个位移。位移是离散的数据，也可以拟合出函数。
弹簧模型可以在sixDof里实现，通过定义pointDisplacement文件中的边界类型为如下

type            sixDoFRigidBodyDisplacement;

可以实现弹簧的定义。
但是问题在于特定的时间范围上指定一个位移，我想到的方法也是在pointDisplacement文件中自定义边界条件实现

type            codedFixedValue;

这样分开定义都可以实现，叠加在一起就矛盾了。各位老师有相关的经验吗，是否有其他的实现方案。

@李东岳
基本理论是通过增加时间步长来尝试滤波，大致的时间步长对应的CFL在5左右，所以我想先用pimple试一下。
这里不是CFL=0.01，是实际应用时CFL只能在1以下，再提高CFL就会发散了。这个0.01是CFL=1对应的时间步长。另外想请问老师设置较大的CFL有什么要求吗，我是直接在controlDict里设置了maxCo。

@李东岳
谢谢老师，那是我的计算设置有问题吗。
文章链接贴在下面
http://dx.doi.org/10.4236/ojfd.2015.52021

各位老师好，

文献中提到两种求解不稳定流动稳态的方法，一个是直接求解稳态方程；另一个是使用瞬态求解器求解，但是增大时间步长达到滤波的效果，也能达到求解稳态的目的。

我在计算不稳定流动的稳态解时尝试使用直接求解稳态的方案，可以较好的计算出对称结构的稳态解（例如simple求解Re=60的圆柱绕流）。但是使用瞬态求解时，pimple方法仍然受限于CFL数，很难增大时间步长。
我的网格最小尺寸大概在0.01，入口速度1，CFL限制下的deltaT大约在0.01

实际操作下确实最大时间步长不能超过0.008，过大就马上发散了。
但是pimple方法不是可以有很大的时间步长和库朗数吗，是不是我设置有问题，controlDict贴在下面。deltaT 0.05;是尝试达到的时间步长，但是是发散的。

application     pimpleFoam;

startFrom       latestTime;

startTime       0;

stopAt          endTime;

endTime         150;

deltaT          0.05;

writeControl    timeStep;

writeInterval   1000;

purgeWrite      50;

writeFormat     binary;

writePrecision  8;

writeCompression off;

timeFormat      general;

timePrecision   8;

runTimeModifiable true;

adjustTimeStep  yes;

maxCo           20;

还有pimple的设置

PIMPLE
{
    correctPhi          yes;
    nOuterCorrectors    50;
    nCorrectors         2;
    nNonOrthogonalCorrectors 0;
	pRefCell 1001; 
    pRefValue 0; 
}
relaxationFactors
{
    fields
    {
    }
    equations
    {
        "U.*"           1;
    }
}

@李东岳
谢谢老师，一开始我认为是入口的速度梯度计算的不准确，所以用差分方法直接计算。后面对照下看来是接近的。
我认为可能是不可压缩的算法里不考虑严格的能量守恒，所以 直接计算并不为0.
但是上面公式推导中使用高斯方法转换体面积分的时候默认了不可压缩时这项为零。应该写出完整的形式计算

Total viscous power in the domain (volume integral): -1.50451907
Total integrated tau : gradU: 1.48509576
Total integrated viscous power (surface): 0.00384099619

这样计算后数值是接近的，等稳定后误差应该还会减小。

@李东岳
老师好，我分开回答
1.这里略去的原因是
在总的能量方程中：

如果不考虑热源、辐射、外力和温度变化。方程可以写成

再考虑动量方程，就得到

因为不可压缩，所以这项就为零，省略去了。
2.离散的形式我没有自己展开，使用的下面代码，这是体积分

// 计算应变速率张量的偏量（deviatoric part of the strain rate tensor）
volTensorField devD = dev(gradU);

// 计算粘性应力张量
volTensorField tau = mu * devD;

// 计算粘性功率密度
volScalarField viscousPowerDensity = U & fvc::div(tau);

// 计算整个计算域的粘性功率（体积分）
scalar totalViscousPowerVolume = 0.0;
forAll(viscousPowerDensity, cellI)
{
    scalar cellVolume = mesh.V()[cellI];
    totalViscousPowerVolume += viscousPowerDensity[cellI] * cellVolume;
}

Info << "Total viscous power in the domain (volume integral): " << totalViscousPowerVolume << endl;

这部分是面积分，有点长

// 计算粘性功率（面积分）
scalar totalViscousPowerSurface = 0.0;

// 定义入口ID
label inletPatchID = mesh.boundaryMesh().findPatchID("inlet"); // 假设入口补丁名为 "inlet"

forAll(mesh.boundaryMesh(), patchI)
{
    const polyPatch& patch = mesh.boundaryMesh()[patchI];
    const word patchName = patch.name();

    const fvPatchVectorField& U_patch = U.boundaryField()[patchI];
    const fvPatchTensorField& tau_patch = tau.boundaryField()[patchI];

    scalar patchViscousPower = 0.0;

    forAll(patch, faceI)
    {
        const vector& Sf = patch.faceAreas()[faceI];
        const vector& U_face = U_patch[faceI];
        tensor tau_face;

        if (patchI == inletPatchID)
        {
            // 计算入口速度梯度的估计
            label cellI = patch.faceCells()[faceI];

            // 使用有限差分法计算速度梯度
            scalar dx = mesh.V()[cellI] / mag(Sf); // 单元在面法向量方向上的长度
            tensor gradU_approx = tensor::zero;
            gradU_approx.xx() = 2 * U_face.x() / dx; // x 方向上的速度梯度
            gradU_approx.yy() = 2 * U_face.y() / dx; // y 方向上的速度梯度
            gradU_approx.zz() = 2 * U_face.z() / dx; // z 方向上的速度梯度
            //Info << "dx: " << dx << endl;
            tau_face = (2.0/3.0) * mu * gradU_approx;
        }
        else
        {
            tau_face = tau_patch[faceI]; // 正常计算
        }

        scalar viscousPowerFace = (tau_face & U_face) & Sf;
        patchViscousPower += viscousPowerFace;

        // 输出壁面上的速度梯度和粘性应力张量的值
        if (patchName == "inlet")
        {
            Info << "patch: " << patchName << ", face: " << faceI << endl;
            Info << "U_face: " << U_face << endl;
            Info << "tau_face: " << tau_face << endl;
            Info << "Sf (face normal vector): " << Sf << endl;
			
            Info << "viscousPowerFace: " << viscousPowerFace << endl;
        }
    }

    totalViscousPowerSurface += patchViscousPower;
    Info << "Patch: " << patchName << ", Integrated viscous power: " << patchViscousPower << endl;
}

@李东岳 抱歉老师，之前的写法有问题，重新补充一下。

这里也是有相同的问题，如果入口使用了充分发展的速度条件，管道流稳定后沿流方向不存在速度梯度，按照面积分的方法求得的结果应该也是0.

各位老师好，

我在计算管道流粘性功率时使用两种积分方案得到的结果误差很大。
在数学上，体积分可以转换为面积分。

我先用整个计算域内的网格计算左边的体积分，得到的值很接近预测值。
再使用右侧的面积分计算几个边界，却得到很小的值。管道流动如下：
左侧入口：充分发展的速度入口，计算得到的面积分很小（接近0）
右侧出口：给定压力出口（0），他的速度梯度是0，所以在这个边界上的面积分是0，和计算得到的值匹配
两侧壁面：无滑移壁面 壁面速度为0，所以边界上的面积分也是0，和计算也相符

问题是出口和壁面都是由于给定的边界条件导致的数值为0，那么入口的面积分结果应该很接近左侧的体积分值才对。

我自己分析是由于入口的速度梯度很难准确计算导致的。所以我使用入口右侧一层的网格速度梯度来代替，但是结果仍然很小。请教大家。

@李东岳 好的我简要描述下问题，就是我按照完整的形式推导了不可压缩流动的能量守恒公式。
但是案例结果并不满足。
我对这一项有疑问
我在程序里输出是有值的，但是之前的回复中提到不可压缩流动中应当为零。
这一项是按照下面这个式子，考虑了速度梯度为零，去掉右边第一项积分来的

我发现这个式子中的三项，积分后的形式都是一致的，感到很困惑。

@李东岳 对的，普遍的能量方程推导。可压缩的能量方程应当也适用不可压缩吧，推导完成后我按照不可压缩的连续性方程做了简化，具体在这个贴子https://www.cfd-china.com/topic/7284/%E4%B8%8D%E5%8F%AF%E5%8E%8B%E7%BC%A9%E6%B5%81%E5%8A%A8%E4%B8%AD%E7%9A%84%E8%83%BD%E9%87%8F%E6%96%B9%E7%A8%8B?_=1721812009357

@李东岳 在 openfoam 不可压缩计算中的能量守恒 中说：

这个方程，如果看积分项，就是左边第4个，为什么存在压力的速度对流？不可压缩不存在这个

老师好，这个公式我又重新检查了，对应的是压力梯度项。

这部分展开后，在不可压缩流动中，右侧第一项确实是0，但是第二项积分应该不为零。同时我也很难理解，这两项积分后通过高斯公式转换为面积分不应该是相同的式子吗，怎么会有不一样的结果呢。

谢谢老师回复，我通过代码输出每个时间步对应的量，发现并没有守恒性。如果理论推导没有问题，可能代码的对应出错了，麻烦帮我看一下：

// 计算速度场的梯度  
    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;
	
	// 计算域内的能量变化率（仅考虑动能）
    volScalarField kineticEnergy = 0.5 * rho * magSqr(U);

    scalar totalKineticEnergyOld = 0.0;
    scalar totalKineticEnergyNew = 0.0;

    forAll(kineticEnergy, cellI)
    {
        totalKineticEnergyOld += kineticEnergy.oldTime()[cellI] * mesh.V()[cellI];
        totalKineticEnergyNew += kineticEnergy[cellI] * mesh.V()[cellI];
    }

    scalar energyChangeRate = (totalKineticEnergyNew - totalKineticEnergyOld) / runTime.deltaTValue();

    Info << "Total Kinetic Energy (old): " << totalKineticEnergyOld << " J" << endl;
    Info << "Total Kinetic Energy (new): " << totalKineticEnergyNew << " J" << endl;
    Info << "Kinetic Energy Change Rate in the domain: " << energyChangeRate << " W" << endl;

        runTime.printExecutionTime(Info);
    }