@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);
    }

各位老师好，

我在处理不可压缩能量守恒时，希望通过分析能量守恒中的各项变化来分析能量流动方向，但是碰到几个问题，请大家指导。

参考《计算流体力学基础及其应用》推导能量方程：

流体微团的能量变化率=流入的净热流+体积力功率+表面力功率

体积力

可以表示为

表面力

分为正压力和切应力，只考虑一个方向的分量

综合三个方向，简化表示为

其中第一项表示可逆的做功过程，第二项为粘性耗散。进一步的这两项各可以拆为两项

净热流

分为热源（辐射等）和传导

能量变化率

这里取守恒形式，将物质导数展开

所以总的能量方程可写为

现在考虑不可压缩流动中，没有外力和热源及温度变化，方程退化为如下形式

将其化简并展开

这里我们比较动量方程和连续性方程

可以消掉能量方程的部分，得到如下表达

也就是说不可压缩流动中，如果不考虑温度变化，粘性耗散项没有作用是能量守恒的基础。

那么当我们考虑这样一个具体问题，不可压缩的圆柱涡激振动，能量守恒的物理意义应该如何表达，是否可以写成如下形式

在全场积分即

我的问题是：
1.这个积分式子的物理意义是
计算域内的能量变化率=所有边界压力功率+所有边界剪切应力功率

还是说

计算域内的能量变化率=所有控制体积上的压力功率+所有控制体积上的剪切应力功率

2.这样就不存在分析各项变化来找到能量流动方向了，主要希望分析流场和结构间的能量流动方向。

@李东岳 谢谢老师，对我帮助很大。
我整理了一下案例和说明放到仓库了https://gitee.com/iceflowww/openfoam-testing.git

@李东岳 那太好了，谢谢老师！

@李东岳 谢谢老师，案例挂在下面了。cavitytest.zip

谢谢老师，我重新调整了案例，麻烦再帮我看一下。
绘制四网格一维模型，并且设置左侧入口速度为1，右侧出口零梯度。

边界条件设置如下

U

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (1 0 0);

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           $internalField;
    
    }
    outlet
    {
        type            zeroGradient;
    }
    frontAndBack
    {
        type            empty;
    }
	Walls
    {
        type            zeroGradient;
    }
}

P

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    inlet
    {
        type            zeroGradient;
    }
    outlet
    {
        type            zeroGradient;
    }
    frontAndBack
    {
        type            empty;
    }
	Walls
    {
        type            zeroGradient;
    }
}

设置对流项的离散为中心差分格式，也叫线性格式

divSchemes
{
    default         none;
    div(phi,U)      Gauss linear;
}

修改动量方程，只保留对流项，并输出各项矩阵系数

 // Momentum predictor
        fvVectorMatrix UEqn
        (
            //fvm::ddt(U)
          fvm::div(phi, U)
          //- fvm::laplacian(nu, U)
        );

        fvVectorMatrix UEqnWithPressure(UEqn == -fvc::grad(p));

        if (piso.momentumPredictor())
        {
            solve(UEqnWithPressure);
        }

        labelUList lAdd = UEqnWithPressure.lduAddr().lowerAddr();
        labelUList uAdd = UEqnWithPressure.lduAddr().upperAddr();
    	labelUList oAdd = UEqnWithPressure.lduAddr().ownerStartAddr();
        scalarField lower = UEqnWithPressure.lower();
        scalarField upper = UEqnWithPressure.upper();
        scalarField diag = UEqnWithPressure.diag();
        vectorField source = UEqnWithPressure.source();

        Info << "lowerAddr:" << lAdd << nl << endl;
        Info << "upperAddr:" << uAdd << nl << endl;
        Info << "ownerStartAddr:" << nl << oAdd << nl << endl;
        Info << "lower:" << lower << nl << endl;
        Info << "upper:" << upper << nl << endl;
        Info << "diag:" << nl << diag << nl << endl;
        Info << "source:" << nl << source << nl << endl;
        //Info << UEqn  << endl;

输出结果如下

Starting time loop

Time = 0.001

Courant Number mean: 0.001 max: 0.001
smoothSolver:  Solving for Ux, Initial residual = 0, Final residual = 0, No Iterations 0
smoothSolver:  Solving for Uy, Initial residual = 0, Final residual = 0, No Iterations 0
lowerAddr:3(0 1 2)

upperAddr:3(1 2 3)

ownerStartAddr:
5(0 1 2 3 3)

lower:3{-0.5}

upper:3{0.5}

diag:
4(0.5 0 0 -0.5)

source:
4{(0 0 0)}

=== Matrix CoeffMat before adding contribution of BCs ===
0.5 0.5 0 0 0
-0.5 0 0.5 0 0
0 -0.5 0 0.5 0
0 0 -0.5 -0.5 0

=== Matrix CoeffMat after adding contribution of BCs ===
0.5 0.5 0 0 1
-0.5 0 0.5 0 0
0 -0.5 0 0.5 0
0 0 -0.5 0.5 0

可以看到，系数矩阵并不是一个对角占优的矩阵，所以无法推进（个人看法）

进一步按照理论计算矩阵内的各项系数：

由于修改后的方程只含有对流项，所以离散结果可以进一步简化：

计算A10如下：

看起来没有问题，但是根据边界条件，只有入口边界的速度为1，如下图

所以 ,按照这样计算的$A_{10}$应该是 -0.25，和输出的矩阵系数也不一样。

进一步的，计算对角矩阵系数（其实和次对角一样），怎么会有0出现，请老师帮忙解惑。

各位老师好，

最近在学习icoFoam，参考李老师写的 CFD:不可压+瞬态算法，有几个问题没有弄清楚，发出来向大家请教。
我按照原生的cavity案例运行，为了计算方便修改成四个网格，并且把边界条件都去掉了。

dimensions      [0 1 -1 0 0 0 0];

internalField   nonuniform List<vector> 4((1 0 0)(1 0 0)(2 0 0)(2 0 0));

boundaryField
{
    movingWall
    {
        type            noSlip;
    }
    fixedWalls
    {
        type            noSlip;
    }
    frontAndBack
    {
        type            empty;
    }
}

主要希望方便计算系数.但是这里就有一个问题，我直接定义了四个网格的速度，我原本以为速度应该处于网格中心，结果看起来处于上下边界。所以直接把所有网格的速度都改为1了，如下图。
根据理论推导的离散结果，可以写出系数矩阵A

这里我输出没有计算之前的矩阵，是一个对角占优的阵，形式应该没问题

=== Matrix CoeffMat before adding contribution of BCs ===
0.0056        -0.0001         -0.0001         0 
-0.0001       0.0056               0          -0.0001 
-0.0001            0                0.0056    -0.0001 
0                  -0.0001          -0.0001   0.0056 

根据上面的公式非对角的元素可以按照位置计算，比如A10（第二行第一列）

完全没对上，不知道哪里有问题，请大家帮忙看下。