原始问题

我老板做模拟，使用 maxCo=1.5 虽然这部分已经是错了，但是神奇的一点是，他的模拟作出来的颗粒的
位置就是不对的。然后我运行 maxCo=0.5 跑出来的颗粒的位置看起来就没有什么问题。

然后我早期的印象中是颗粒会卡在碰到的界面上，然后等下一个Euler时间步。
然后我查了一下源代码，问了一下AI，因为还是没有看懂。。
最后得到的解释是下面的一坨。

OpenFOAM parcel 碰撞处理

User
template
template
void Foam::KinematicParcel::hitFace(TrackData& td)
{
typename TrackData::cloudType::parcelType& p =
static_cast<typename TrackData::cloudType::parcelType&>(*this);

td.cloud().functions().postFace(p, p.face(), td.keepParticle);

}

template
void Foam::KinematicParcel::hitFace(int& td)
{}

template
template
bool Foam::KinematicParcel::hitPatch
(
const polyPatch& pp,
TrackData& td,
const label patchi,
const scalar trackFraction,
const tetIndices& tetIs
)
{
typename TrackData::cloudType::parcelType& p =
static_cast<typename TrackData::cloudType::parcelType&>(*this);

// Invoke post-processing model
td.cloud().functions().postPatch
(
    p,
    pp,
    trackFraction,
    tetIs,
    td.keepParticle
);

// Invoke surface film model
if (td.cloud().surfaceFilm().transferParcel(p, pp, td.keepParticle))
{
    // All interactions done
    return true;
}
else
{
    // Invoke patch interaction model
    return td.cloud().patchInteraction().correct
    (
        p,
        pp,
        td.keepParticle,
        trackFraction,
        tetIs
    );
}

}

template
template
void Foam::KinematicParcel::hitProcessorPatch
(
const processorPolyPatch&,
TrackData& td
)
{
td.switchProcessor = true;
}

template
template
void Foam::KinematicParcel::hitWallPatch
(
const wallPolyPatch& wpp,
TrackData& td,
const tetIndices&
)
{
// Wall interactions handled by generic hitPatch function
}

你看这部分程序，在OpenFOAM中，如果parcel 碰到cell边界会怎么处理？我知道这个程序中有2个时间步，一个是Euler 一个是Lagrangian时间步，如果碰到任意边界的话，颗粒会在哪个时间步中离开当前cell？

Copilot
The user is asking about how parcels (particles) are handled when they hit cell boundaries in OpenFOAM’s Lagrangian tracking framework. They specifically want to know about the two time steps (Euler and Lagrangian) and at which time step the particle leaves the current cell when hitting a boundary.

Let me analyze the code carefully.

Looking at the move function (which I can see in the full file), the tracking logic is:

The move function contains the main tracking loop
Inside the loop, trackToFace is called which moves the parcel to the face
When a faceThe user is asking about how parcels (particles) are handled when they hit cell boundaries in OpenFOAM’s Lagrangian tracking framework. They specifically want to know about the two time steps (Euler and Lagrangian) and at which time step the particle leaves the current cell when hitting a boundary.

Let me analyze the code carefully.

Looking at the move function (which I can see in the full file), the tracking logic is:

The move function contains the main tracking loop
Inside the loop, trackToFace is called which moves the parcel to the face
When a face is hit, hitFace and hitPatch are called
The Eulerian time step is the global CFD time step, while the Lagrangian time step is a sub-step within it, controlled by a Courant number condition for particles. The code calculates the maximum allowable Lagrangian time step based on the Courant number to ensure stability.
In the tracking loop, the code computes the actual time step by taking the minimum of the maximum allowable step and the remaining time, then adjusts it further based on the particle’s velocity and the cell length scale to enforce the Courant condition locally.

Adjusting particle time step… the particle actually leaves its current cell during the Lagrangian sub-step itself, not waiting until the end of the Eulerian time step. After trackToFace completes, if the particle is on a boundary, the patch interaction gets handled through hitPatch or similar functions, which let the surface film model and patch interaction model determine what happens next. The loop then decrements tEnd by the actual time used, so the remaining time continues tracking within the same Eulerian step. is that within each Lagrangian sub-step, the particle moves toward a target position based on its velocity, and if it hits a face before reaching that target, the sub-step is truncated and the particle immediately transitions to the adjacent cell. This process repeats until the full Eulerian time step is consumed, so the particle’s cell change occurs at the precise moment of face crossing during Lagrangian tracking, not at the end of the time step.## OpenFOAM Lagrangian 颗粒碰到边界的处理流程

两个时间步的关系

先看 move 函数中的时间步定义：

scalar tEnd = (1.0 - p.stepFraction())*trackTime;  // 剩余 Euler 时间
scalar dtMax = trackTime;
if (td.cloud().solution().transient())
{
    dtMax *= maxCo;  // Lagrangian 子步上限 = Euler## OpenFOAM Lagrangian 颗粒碰到边界的处理流程

### 两个时间步的关系

先看 `move` 函数中的时间步定义：

```cpp
scalar tEnd = (1.0 - p.stepFraction())*trackTime;  // 剩余 Euler 时间
scalar dtMax = trackTime;
if (td.cloud().solution().transient())
{
    dtMax *= maxCo;  // Lagrangian 子步上限 = Euler步 × Courant数
}

时间步	变量	说明
Euler 时间步	`trackTime`	全局 CFD 时间步长
Lagrangian 子步	`dt`（≤ `dtMax`）	由颗粒 Courant 数限制的子步

颗粒碰到边界的关键时刻：Lagrangian 子步内

颗粒在 Lagrangian 子步内离开当前 cell，具体发生在 trackToFace 调用中：

// 在 while 循环的每个 Lagrangian 子步中：
const scalar d = dt*magU;
const scalar dCorr = min(d, maxCo*cellLengthScale[celli]);
dt *=
    dCorr/d
   *p.trackToFace(p.position() + dCorr*U_/magU, td);  // ← 颗粒在此穿越面

trackToFace 返回一个 0~1 的分数（实际走完的比例）。若中途命中面，则：

返回值 < 1，dt 被按比例缩短
内部自动触发 hitFace() → hitPatch()

碰到边界后的处理链（`hitPatch` 中）

触发 postPatch()（后处理函数）
        ↓
surfaceFilm().transferParcel()   ← 先判断是否被液膜吸收
    是 → return true（颗粒消失，keepParticle 由 film 控制）
    否 ↓
patchInteraction().correct()     ← 再判断壁面交互（反弹/吸附/穿透）

hitWallPatch 本身为空，说明壁面交互完全由 hitPatch 中的 patchInteraction 统一负责。

结论

颗粒在 Lagrangian 子步内（trackToFace 调用时） 离开当前 cell，并非等到 Euler 时间步结束。之后：

若 keepParticle = false（被吸附/消失）→ while 循环退出
若 switchProcessor = true（跨处理器）→ while 循环退出
若反弹/穿越内部面 → 更新 cell，tEnd 减去已用 dt，继续在同一 Euler 时间步内的 while 循环追踪剩余时间

User
什么意思，你在给我详细的解释一下，并根据碰到了不同的位置，而有不同的操作解释。如果发生跨face和处理器是分别怎么处理的？

Copilot
The user wants a more detailed explanation of what happens when a parcel hits different types of boundaries in OpenFOAM. They specifically want to understand:

What happens when crossing an internal face (between cells)
What happens when crossing a processor patch (parallel decomposition boundary)
What happens when hitting a wall patch

Let me trace through the code carefully.

The key function is trackToFace (called in move), which internally calls hitFace and hitPatch. Let me look at the particle.C base class to understand trackToFace better, but based on the code provided, I can explain the flow.
When a particle crosses an internal face between cells, hitFace gets called to handle post-processing, but hitPatch is skipped since it’s not a boundary patch—the particle just moves into the adjacent cell and the tracking loop continues. At a processor patch boundary, both hitFace and hitPatch are invoked, which is where the particle gets transferred to the appropriate processor for parallel computation. The while loop exits when td.switchProcessor becomes true, leaving the particle positioned at the processor patch boundary where the cloud’s parallel communication takes over the transfer. There’s also an explicit check after tracking that sets td.switchProcessor if the particle lands on a processor patch, though hitProcessorPatch may already handle this—it could be redundant or cover an edge case. Now looking at what happens when a particle hits a wall patch instead. For other patches without special handling, the particle is simply removed by setting keepParticle = false. Now I’m tracing through the actual call sequence in OpenFOAM’s particle tracking: trackToFace advances the particle, then when a face is crossed it triggers hitFace, which checks if it’s a boundary face and calls the appropriate hitPatch method based on the patch type. Now I’ll walk through each scenario—internal face crossing, wall patches, processor patches, and other boundaries—to show how the particle tracking handles them differently.## 颗粒碰到不同位置的详细处理流程

先建立整体调用链的认知：

move() 的 while 循环
    └─ trackToFace()          ← 几何追踪，遇到面就停
           └─ hitFace()       ← 每次穿越任意面都调用
                  └─ hitPatch() / hitWallPatch() / hitProcessorPatch()
                               ← 只在边界面时调用（内部面不调用）

情形一：穿越内部面（cell → 相邻 cell）

// trackToFace 内部发现穿越了内部面，调用：
hitFace(TrackData& td)
{
    // 仅触发后处理回调，无其他操作
    td.cloud().functions().postFace(p, p.face(), td.keepParticle);
}

没有 hitPatch 调用。trackToFace 返回后，dt 被缩短为实际走到面的时间比例，tEnd 减去该 dt，while 循环继续，颗粒在新 cell 中完成剩余追踪。

1 2	while循环第N次: [旧cell] ──trackToFace──> 碰到内部面 → hitFace → dt缩短 while循环第N+1次: [新cell] 继续追踪剩余 tEnd

情形二：碰到处理器边界（并行分解面）

调用链：

1
2
3

hitFace()           → postFace 回调
hitPatch(5个参数)   → postPatch + surfaceFilm + patchInteraction
hitProcessorPatch() → td.switchProcessor = true   ← 关键！

然后 move() 中还有一个显式检查：

if (p.onBoundary() && td.keepParticle)
{
    if (isA<processorPolyPatch>(pbMesh[p.patch(p.face())]))
    {
        td.switchProcessor = true;  // 双重保险
    }
}

while 循环的退出条件是 !td.switchProcessor，一旦为 true，循环立即退出。颗粒停在处理器边界面上，由 Cloud 层的并行通信（KinematicCloud::motion() 中的 transfer()）将其序列化发送到邻居处理器，在那边以剩余的 stepFraction 继续追踪。

处理器A: while循环 → 碰到processorPatch → switchProcessor=true → 退出循环
                                                                       ↓
                                                         Cloud::transfer() 发送颗粒
                                                                       ↓
处理器B: move() 重新启动, tEnd = (1 - stepFraction)*trackTime 继续追踪

情形三：碰到壁面（wallPolyPatch）

调用链：

1
2
3

hitFace()           → postFace 回调
hitPatch(5个参数)   → ① postPatch ② surfaceFilm ③ patchInteraction
hitWallPatch()      → 什么都不做（注释已说明）

核心在 hitPatch 中：

// ① 后处理回调
td.cloud().functions().postPatch(p, pp, trackFraction, tetIs, td.keepParticle);

// ② 液膜模型优先判断
if (td.cloud().surfaceFilm().transferParcel(p, pp, td.keepParticle))
{
    return true;  // 颗粒被液膜吸收，结束
}
else
{
    // ③ 壁面交互模型决定命运
    return td.cloud().patchInteraction().correct(p, pp, td.keepParticle, ...);
}

patchInteraction 模型（如 rebound、standardWallInteraction）会修改速度并设置 keepParticle：

patchInteraction 结果	keepParticle	while 循环
反弹（rebound）	`true`	继续，以反弹后速度追踪
吸附/沉积（stick）	`false`	退出循环，颗粒死亡
逃逸（escape）	`false`	退出循环，颗粒死亡

情形四：碰到其他边界面（非壁面、非处理器）

触发 2 参数版本的 fallback hitPatch：

template<class ParcelType>
template<class TrackData>
void Foam::KinematicParcel<ParcelType>::hitPatch
(
    const polyPatch&,   // 任意未特殊处理的 patch
    TrackData& td
)
{
    td.keepParticle = false;  // 直接杀死颗粒
}

while 循环因 !td.keepParticle 退出，颗粒被移除。

汇总对比

颗粒碰到不同位置的详细处理流程

先建立整体调用链的认知：

move() 的 while 循环
    └─ trackToFace()          ← 几何追踪，遇到面就停
           └─ hitFace()       ← 每次穿越任意面都调用
                  └─ hitPatch() / hitWallPatch() / hitProcessorPatch()
                               ← 只在边界面时调用（内部面不调用）

情形一：穿越内部面（cell → 相邻 cell）

// trackToFace 内部发现穿越了内部面，调用：
hitFace(TrackData& td)
{
    // 仅触发后处理回调，无其他操作
    td.cloud().functions().postFace(p, p.face(), td.keepParticle);
}

没有 hitPatch 调用。trackToFace 返回后，dt 被缩短为实际走到面的时间比例，tEnd 减去该 dt，while 循环继续，颗粒在新 cell 中完成剩余追踪。

1 2	while循环第N次: [旧cell] ──trackToFace──> 碰到内部面 → hitFace → dt缩短 while循环第N+1次: [新cell] 继续追踪剩余 tEnd

情形二：碰到处理器边界（并行分解面）

调用链：

1
2
3

hitFace()           → postFace 回调
hitPatch(5个参数)   → postPatch + surfaceFilm + patchInteraction
hitProcessorPatch() → td.switchProcessor = true   ← 关键！

然后 move() 中还有一个显式检查：

if (p.onBoundary() && td.keepParticle)
{
    if (isA<processorPolyPatch>(pbMesh[p.patch(p.face())]))
    {
        td.switchProcessor = true;  // 双重保险
    }
}

处理器A: while循环 → 碰到processorPatch → switchProcessor=true → 退出循环
                                                                       ↓
                                                         Cloud::transfer() 发送颗粒
                                                                       ↓
处理器B: move() 重新启动, tEnd = (1 - stepFraction)*trackTime 继续追踪

情形三：碰到壁面（wallPolyPatch）

调用链：

1
2
3

hitFace()           → postFace 回调
hitPatch(5个参数)   → ① postPatch ② surfaceFilm ③ patchInteraction
hitWallPatch()      → 什么都不做（注释已说明）

核心在 hitPatch 中：

// ① 后处理回调
td.cloud().functions().postPatch(p, pp, trackFraction, tetIs, td.keepParticle);

// ② 液膜模型优先判断
if (td.cloud().surfaceFilm().transferParcel(p, pp, td.keepParticle))
{
    return true;  // 颗粒被液膜吸收，结束
}
else
{
    // ③ 壁面交互模型决定命运
    return td.cloud().patchInteraction().correct(p, pp, td.keepParticle, ...);
}

patchInteraction 模型（如 rebound、standardWallInteraction）会修改速度并设置 keepParticle：

patchInteraction 结果	keepParticle	while 循环
反弹（rebound）	`true`	继续，以反弹后速度追踪
吸附/沉积（stick）	`false`	退出循环，颗粒死亡
逃逸（escape）	`false`	退出循环，颗粒死亡

情形四：碰到其他边界面（非壁面、非处理器）

触发 2 参数版本的 fallback hitPatch：

template<class ParcelType>
template<class TrackData>
void Foam::KinematicParcel<ParcelType>::hitPatch
(
    const polyPatch&,   // 任意未特殊处理的 patch
    TrackData& td
)
{
    td.keepParticle = false;  // 直接杀死颗粒
}

while 循环因 !td.keepParticle 退出，颗粒被移除。

汇总对比

碰到内部面         → hitFace → while继续，进入新cell
碰到processorPatch → hitFace → hitPatch → hitProcessorPatch → switchProcessor=true → 跨处理器传输
碰到wallPatch      → hitFace → hitPatch → [液膜? / patchInteraction决定反弹or死亡]
碰到其他patch      → hitFace → hitPatch(fallback) → keepParticle=false → 颗粒消失

User
那么如果遇到跨处理器传输的话，涉及到Euler时间步么？

Copilot
The user is asking whether the Euler time step is involved when a particle crosses a processor boundary (parallel decomposition). Let me think about this carefully.

When a particle hits a processor patch:

hitProcessorPatch sets td.switchProcessor = true
The while loop in move() exits because of !td.switchProcessor
The particle is at the processor boundary face with a certain stepFraction

The key question is: when the particle is transferred to the neighboring processor and resumes tracking, does it use the same Euler time step or a new one?

Looking at the code:

1	scalar tEnd = (1.0 - p.stepFraction())*trackTime;

stepFraction records how much of the Euler time step has been consumed. When the particle is transferred to the other processor, it carries its stepFraction with it. On the receiving processor, move() is called again with the same trackTime (Euler time step), and tEnd is calculated from the remaining stepFraction. So the Euler time step is NOT reset—the particle continues within the same Euler time step on the receiving processor, using only the remaining time. This is exactly what stepFraction is designed for: allowing particles to be transferred mid-Euler-step while maintaining continuity.

The mechanism works by updating stepFraction after each Lagrangian sub-step to reflect progress through the Euler step, then when a particle crosses a processor boundary, it carries this updated value to the neighboring processor where move() resumes with the remaining time calculated from the fraction.

1
2
3


cpp
scalar tEnd = (1.0 - p.stepFraction())*trackTime;

This ensures the particle continues tracking through the same Euler time step seamlessly across processor boundaries, maintaining consistency in the Lagrangian tracking throughout the parallel computation.## 跨处理器传输与 Euler 时间步的关系

不会重新开始新的 Euler 时间步。 关键在于 stepFraction。

stepFraction 的作用

在 move() 中，每个 Lagrangian 子步结束后都会更新 stepFraction：

scalar tEnd = (1.0 - p.stepFraction())*trackTime;  // 用已记录进度算剩余时间
// ...（追踪一个子步）...
scalar newStepFraction = 1.0 - tEnd/trackTime;      // 更新进度
p.stepFraction() = newStepFraction;                 // 存入颗粒本身

stepFraction 记录的是当前 Euler 步已完成的比例（0→1），它跟随颗粒数据一起被序列化传输到邻居处理器。

传输后的恢复

邻居处理器收到颗粒后，重新调用 move(trackTime) ——注意 trackTime 是同一个全局 Euler 时间步长：

// 邻居处理器上，用传过来的 stepFraction 计算剩余时间
scalar tEnd = (1.0 - p.stepFraction()) * trackTime;
//             └──────────────────────┘
//             不是 0，而是从中断处继续

如果传输前 stepFraction = 0.6，那么邻居处理器只需追踪剩余的 40% 时间，而不是从头开始。

完整的跨处理器时序

Euler时间步 T ~ T+dt（全局统一）
│
├── 处理器A: move(dt)
│     ├── 子步1: stepFraction 0.0 → 0.3  [正常追踪]
│     ├── 子步2: stepFraction 0.3 → 0.6  [正常追踪]
│     └── 子步3: 碰到processorPatch
│                 switchProcessor=true → 退出while循环
│                 stepFraction 已更新为 0.6
│                 颗粒（携带stepFraction=0.6）→ 传输 →
│
└── 处理器B: move(dt)   ← 同一个dt！
      tEnd = (1-0.6)*dt = 0.4*dt   ← 从断点继续
      ├── 子步4: stepFraction 0.6 → 0.85  [继续追踪]
      └── 子步5: stepFraction 0.85 → 1.0  [完成整个Euler步]

结论

跨处理器传输是在 Lagrangian 子步级别中断并恢复的，Euler 时间步作为整体跨越了处理器边界，颗粒通过 stepFraction 保持在同一 Euler 步内的连续性，传输本身对 Euler 时间步没有任何影响。

User
但是有2个模拟，他们的Euler 时间步是和maxCo相关的，如果maxCo比较大 1.5 和maxCo0.5的到的结果就完全不一样，为什么会这样？同样的模拟，只有maxCo不一样

Copilot
The user is asking why two simulations with the same setup but different maxCo values (1.5 vs 0.5) give completely different results. This is a good question about numerical accuracy in OpenFOAM’s Lagrangian tracking.

Let me analyze the code carefully to understand how maxCo affects the simulation.

Looking at the move() function:

scalar dtMax = trackTime;
if (td.cloud().solution().transient())
{
    dtMax *= maxCo;  // Lagrangian sub-step upper limit
}

And in the tracking loop:

1	const scalar d = dt*magU;

The code limits particle displacement per sub-step by computing the corrected distance as the minimum of the actual distance traveled and the maximum allowed distance based on the cell length scale and maxCo. This maxCo value acts as the particle Courant number, controlling both the maximum sub-step duration and the maximum displacement within each cell to ensure numerical stability.

The reason different maxCo values produce different results comes down to how forces are integrated during each sub-step. When calculating particle velocity, the integration scheme uses the Lagrangian sub-step size dt to advance the particle’s equation of motion. A larger maxCo means larger sub-steps, which changes how the force integration accumulates over time and ultimately affects the particle trajectory.

The critical issue is that particles move based on their velocity from the previous sub-step, then only after tracking to the new position do we update the fluid properties and recalculate the velocity for the next step. This lag between movement and force recalculation means that larger sub-steps can accumulate more error, especially in non-uniform flow fields.

The dCorr correction factor limits how far a particle can travel in a single sub-step based on the maxCo number and cell length scale, so when maxCo = 1.5, the correction becomes more aggressive in constraining the step size.

With larger maxCo values, particles can traverse multiple cells in one sub-step without sampling the fluid properties in between, leading to incorrect drag calculations and potentially missing interaction regions entirely. This is fundamentally a numerical stability problem with the explicit integration scheme — the larger time steps make the force integration less accurate. the old cell index is still passed to calc(), which might use stale cell-based properties. The interpolation in setCellValues does use the current position after crossing, but there’s a mismatch between what position-based values are computed and what cell-based values calc() receives. particle has moved far away — the momentum transfer goes to the wrong cell. There are also boundary interactions to consider: with large maxCo, particles might reach boundaries in fewer sub-steps with different accumulated forces. The core insight is that maxCo controls the accuracy of Lagrangian tracking — smaller values give more accurate results but cost more computationally, and the results should converge as maxCo decreases rather than remain identical.## 为什么 maxCo 不同会导致结果完全不同

这是一个典型的数值精度问题，根源在于 move() 的执行顺序有一个关键的”先移动，再算力”的结构缺陷。

核心问题：位置更新与速度更新的错位

仔细看 while 循环里的执行顺序：

const label celli = p.cell();        // ① 记录当前cell（旧位置）

// ② 用【旧速度 U_】做几何追踪，移动颗粒到新位置
dt *= dCorr/d * p.trackToFace(p.position() + dCorr*U_/magU, td);

tEnd -= dt;
// ...

// ③ 用【新位置】的流场属性更新 Uc_, rhoc_...
p.setCellValues(td, dt, celli);

// ④ 用新流场属性计算【新速度 U_】，并把动量源项累加到【旧 celli】
p.calc(td, dt, celli);

每个 Lagrangian 子步都是显式的：颗粒用旧速度走一步，然后才感受到新位置的流场并更新速度。

maxCo 如何放大这个误差

maxCo 控制每个子步最多走多远：

1	const scalar dCorr = min(d, maxCo * cellLengthScale[celli]);

maxCo	每子步最大位移	子步数	累积误差
0.5	0.5 个 cell 长度	多	小
1.5	1.5 个 cell 长度	少	大

当 maxCo = 1.5 时，颗粒一步可以跨越超过 1 个 cell，带来 3 个具体问题：

问题一：跨 cell 采样缺失

1 2	maxCo=0.5: A → B → C → D (4步，每个cell都采样了流场) maxCo=1.5: A → C → D (2步，B cell的流场被完全跳过)

颗粒在 B cell 中没有感受到任何力的变化，拖曳力计算出错。

问题二：动量源项累积到错误的 cell

calc() 把双向耦合的动量源项累积到 celli（子步开始时的 cell）：

1	td.cloud().UTrans()[celli] += np0*dUTrans; // celli 是旧cell！

当 maxCo=1.5 时颗粒已经跑到 1.5 个 cell 之外，动量交换仍算在起点 cell，导致流场反馈位置错误。

问题三：力积分的显式误差

calcVelocity 中对运动方程进行时间积分，步长是 dt：

$$U_{new} = \text{integrate}(U_{old},\ dt,\ \underbrace{a(x_{new})}_{\text{新位置的加速度}})$$

为什么 maxCo 不同会导致结果完全不同

这是一个典型的数值精度问题，根源在于 move() 的执行顺序有一个关键的”先移动，再算力”的结构缺陷。

核心问题：位置更新与速度更新的错位

仔细看 while 循环里的执行顺序：

const label celli = p.cell();        // ① 记录当前cell（旧位置）

// ② 用【旧速度 U_】做几何追踪，移动颗粒到新位置
dt *= dCorr/d * p.trackToFace(p.position() + dCorr*U_/magU, td);

tEnd -= dt;
// ...

// ③ 用【新位置】的流场属性更新 Uc_, rhoc_...
p.setCellValues(td, dt, celli);

// ④ 用新流场属性计算【新速度 U_】，并把动量源项累加到【旧 celli】
p.calc(td, dt, celli);

每个 Lagrangian 子步都是显式的：颗粒用旧速度走一步，然后才感受到新位置的流场并更新速度。

maxCo 如何放大这个误差

maxCo 控制每个子步最多走多远：

1	const scalar dCorr = min(d, maxCo * cellLengthScale[celli]);

maxCo	每子步最大位移	子步数	累积误差
0.5	0.5 个 cell 长度	多	小
1.5	1.5 个 cell 长度	少	大

当 maxCo = 1.5 时，颗粒一步可以跨越超过 1 个 cell，带来 3 个具体问题：

问题一：跨 cell 采样缺失

1 2	maxCo=0.5: A → B → C → D (4步，每个cell都采样了流场) maxCo=1.5: A → C → D (2步，B cell的流场被完全跳过)

颗粒在 B cell 中没有感受到任何力的变化，拖曳力计算出错。

问题二：动量源项累积到错误的 cell

calc() 把双向耦合的动量源项累积到 celli（子步开始时的 cell）：

1	td.cloud().UTrans()[celli] += np0*dUTrans; // celli 是旧cell！

当 maxCo=1.5 时颗粒已经跑到 1.5 个 cell 之外，动量交换仍算在起点 cell，导致流场反馈位置错误。

问题三：力积分的显式误差

calcVelocity 中对运动方程进行时间积分，步长是 dt：

$$U_{new} = \text{integrate}(U_{old},\ dt,\ \underbrace{a(x_{new})}_{\text{新位置的加速度}})$$

maxCo 越大 → dt 越大 → 显式积分误差越大，在高曳力区域（如颗粒与流体速度差大时）尤其明显。

结论

maxCo 在 Lagrangian 追踪里不只是 CFL 稳定性条件，更是精度参数。两个模拟给出不同结果是预期内的——它们本质上在用不同精度求解同一个问题：

1
2
3

maxCo → 0  ：理论上的"精确"解
maxCo = 0.5：较精确，结果趋于收敛
maxCo = 1.5：误差较大，结果不可靠

建议做 maxCo 敏感性分析（0.5 → 0.25 → 0.1），当结果不再变化时，该 maxCo 才是足够精确的。对于强曳力或复杂流场，通常需要 maxCo ≤ 0.5。

User
我说的maxCo 不是lagrangian的maxCo，我指的是流场的最大Co，在OpenFoam的controlDict中的那个maxCo 用于控制Euler time step的那个

Copilot
The user is clarifying that they’re talking about the fluid (Eulerian) maxCo in controlDict, which controls the Euler time step size via the adjustTimeStep feature, not the Lagrangian particle Courant number. This is a different parameter entirely.

So the question is: why do two simulations with the same setup but different Eulerian maxCo (controlling the fluid time step size) give different Lagrangian particle results?

Let me think about this carefully.

In OpenFOAM, when adjustTimeStep yes and maxCo is set in controlDict, the solver adjusts deltaT (the Euler time step) to maintain the fluid Courant number below maxCo.

So with a higher maxCo value, the Euler time steps become larger, which means the particleThe user is clarifying that they’re talking about the fluid (Eulerian) maxCo in controlDict, which controls the Euler time step size via the adjustTimeStep feature, not the Lagrangian particle Courant number. This is a different parameter entirely.

So the question is: why do two simulations with the same setup but different Eulerian maxCo (controlling the fluid time step size) give different Lagrangian particle results?

Let me think about this carefully.

In OpenFOAM, when adjustTimeStep yes and maxCo is set in controlDict, the solver adjusts deltaT (the Euler time step) to maintain the fluid Courant number below maxCo.

So with a higher maxCo value, the Euler time steps become larger, which means the particle tracking time step also increases proportionally. Inside the Lagrangian particle tracking code, there’s a separate maxCo parameter that further subdivides the Euler time step for particle integration, but the base tracking time is directly tied to the fluid’s deltaT. The realThe user is clarifying that they’re talking about the fluid (Eulerian) maxCo in controlDict, which controls the Euler time step size via the adjustTimeStep feature, not the Lagrangian particle Courant number. This is a different parameter entirely.

So the question is: why do two simulations with the same setup but different Eulerian maxCo (controlling the fluid time step size) give different Lagrangian particle results?

Let me think about this carefully.

In OpenFOAM, when adjustTimeStep yes and maxCo is set in controlDict, the solver adjusts deltaT (the Euler time step) to maintain the fluid Courant number below maxCo.

So with a higher maxCo value, the Euler time steps become larger, which means the particle tracking time step also increases proportionally. Inside the Lagrangian particle tracking code, there’s a separate maxCo parameter that further subdivides the Euler time step for particle integration, but the base tracking time is directly tied to the fluid’s deltaT. The real question is whether running both simulations to the same physical time produces different results. With a smaller Eulerian Courant number, I get more time steps but smaller individual steps, while a larger Courant number means fewer, larger steps. The Lagrangian particles track through each Euler step with sub-steps scaled to the Euler time step size, so the particle trajectories depend on how the fluid velocity field evolves—and that field only updates at each Euler step, staying frozen in between.

This creates an operator splitting error: with larger Euler time steps, the fluid field used for Lagrangian tracking becomes increasingly stale, so particles are following a less-updated velocity field over longer intervals. The two simulations will diverge because they’re essentially solving the coupled system with different temporal resolutions and different amounts of lag between when the fluid updates and when particles respond to it. -steps magnifies the drag force error since Reynolds number and acceleration are computed from stale carrier velocity data. The integration scheme itself introduces additional error—analytical schemes assume constant acceleration over the sub-step, while Euler schemes have truncation error proportional to dt^2, so larger Euler time steps compound these inaccuracies.

Beyond the sub-step mechanics, there’s operator splitting error where the fluid field remains frozen for longer periods, causing the Lagrangian particles to see a less accurate snapshot of the instantaneous flow. There’s also accumulated source term error from how momentum and heat sources are applied to the fluid in discrete lumps—larger time steps mean larger lumping errors. The core problem is that the Lagrangian-Eulerian coupling itself is only first-order accurate in time, and the Eulerian maxCo directly controls deltaT, which sets the scale of all these coupling errors.## 流场 maxCo（Euler 时间步）影响 Lagrangian 结果的原因

这是一个更根本的问题。Euler deltaT 直接作为 trackTime 传入 Lagrangian 的 move()：

1 2	// 在 KinematicCloud 中调用： p.move(td, deltaT); // deltaT 就是 Euler 时间步

然后在 move() 里：

1
2
3

scalar tEnd = (1.0 - p.stepFraction()) * trackTime;  // trackTime = Euler deltaT
scalar dtMax = trackTime;
if (transient) dtMax *= lagrangianMaxCo;              // Lagrangian 子步上限 ∝ deltaT

Euler deltaT 越大，Lagrangian 子步也等比例变大。 这导致了两层误差来源：

第一层：算子分裂误差（Operator Splitting）

OpenFOAM 的耦合求解是交错进行的：

每个 Euler 步内：
① 流场求解器推进 U_fluid (t → t + deltaT)
② Lagrangian 追踪颗粒，但流场 U_fluid 是冻结的（固定为 t 时刻的值）
③ 颗粒的动量源项累积后一次性施加给流场

这就是算子分裂。流场对颗粒的影响 和 颗粒对流场的反馈 在每个 Euler 步内是解耦的：

$$\underbrace{U_c^{frozen}}{\text{整个 Lagrangian 追踪用这一个值}} \neq \underbrace{U_c(t’)}{\text{实际在 } t’ \in [t, t+\Delta t] \text{ 的瞬时值}}$$

Euler maxCo	deltaT	流场冻结时长	算子分裂误差
0.5	小	短	小
1.5	大	长	大

在非定常流场（如湍流）中，deltaT 越大，颗粒感受到的流场越”过时”，拖曳力计算越不准确。

第二层：动量源项的时间积分误差

双向耦合时，颗粒对流场的动量源项是在整个 Euler 步内累积后一次性施加的：

// 在整个 Euler 步的所有 Lagrangian 子步中累积：
td.cloud().UTrans()[celli] += np0 * dUTrans;

// Euler 步结束后，一次性施加给流场方程（在 UEqn 中）：
// -fvm::Sp(UCoeff/V, U) + UTrans/V

这相当于把一段时间内连续的反馈力集中到一个时刻施加。deltaT 越大，这个近似越粗糙：

1 2	maxCo=0.5 (小deltaT): 每隔很短时间校正一次流场 → 紧密耦合 maxCo=1.5 (大deltaT): 很久才校正一次 → 耦合松散，震荡或发散

第三层：Re 和拖曳力的计算误差

calc() 在每个 Lagrangian 子步开始时计算 Reynolds 数和拖曳力：

1
2
3

const scalar Re = this->Re(U_, d_, rhoc_, muc_);
// 基于此时的 U_ 和 Uc_（已冻结）计算出 Feff
// 然后用 Feff 积分整个子步 dt

子步 dt ∝ deltaT。若颗粒在大 dt 内速度变化显著（如高曳力区），用步初的 Re 和 Feff 来代表整个步长是不准确的。

本质总结

这是一个 Euler-Lagrange 耦合的时间精度问题：

$$\text{总误差} \sim O(\Delta t) \quad \text{（一阶精度的算子分裂）}$$

两个模拟在物理上是同一个问题，但在数值上求解的是不同精度的近似。deltaT 越小，Euler 和 Lagrange 的耦合越紧密，结果越收敛于真实解。

实际上，Lagrangian 模拟对 Euler 时间步的要求往往比纯流场模拟更严格，因为颗粒的响应时间（Stokes number）可能远小于流体的对流时间尺度。

OF4中颗粒碰撞内部面的处理

原始问题

OpenFOAM parcel 碰撞处理

两个时间步的关系

颗粒碰到边界的关键时刻：Lagrangian 子步内

碰到边界后的处理链（hitPatch 中）

结论

情形一：穿越内部面（cell → 相邻 cell）

情形二：碰到处理器边界（并行分解面）

情形三：碰到壁面（wallPolyPatch）

情形四：碰到其他边界面（非壁面、非处理器）

汇总对比

颗粒碰到不同位置的详细处理流程

情形一：穿越内部面（cell → 相邻 cell）

情形二：碰到处理器边界（并行分解面）

情形三：碰到壁面（wallPolyPatch）

情形四：碰到其他边界面（非壁面、非处理器）

汇总对比

stepFraction 的作用

传输后的恢复

完整的跨处理器时序

结论

核心问题：位置更新与速度更新的错位

maxCo 如何放大这个误差

问题一：跨 cell 采样缺失

问题二：动量源项累积到错误的 cell

问题三：力积分的显式误差

为什么 maxCo 不同会导致结果完全不同

核心问题：位置更新与速度更新的错位

maxCo 如何放大这个误差

问题一：跨 cell 采样缺失

问题二：动量源项累积到错误的 cell

问题三：力积分的显式误差

结论

第一层：算子分裂误差（Operator Splitting）

第二层：动量源项的时间积分误差

第三层：Re 和拖曳力的计算误差

本质总结

碰到边界后的处理链（`hitPatch` 中）