SPH结合N体引力模拟中行星旋转阻尼问题及优化问询

问题核心 #

我实现了结合引力与SPH的N体模拟，可模拟行星形变、碰撞等效果，但现有粘性力计算会无差别阻尼刚性旋转的行星（即使粒子间无相对刚体运动）。当前粘性力计算代码：

visc_force = viscosity * Mass_neighbor * (Vel_neighbor - Vel_particle) / density_neighbor * kernelFun_laplacianf(r, r0);

粘性是行星聚集必需的能量耗散机制（纯引力无法形成稳定行星），但刚体旋转不应被阻尼；此前尝试仅保留粘性力径向分量的方案，会破坏层流等真实流体效果。现寻求无需额外循环或扩大邻域的高效真实粘性计算方案。

原有公式的问题分析 #

你使用的是SPH中最简单的拉普拉斯型粘性项，它会对所有惯性系下的速度差施加阻尼——而刚体旋转时，相邻粒子在惯性系下存在切向速度差（尽管随体坐标系下无相对运动），这个项会误将该速度差判定为需要耗散的剪切变形，从而阻尼旋转。

高效解决方案：保留剪切粘性，排除刚体旋转分量 #

核心思路是：仅对粒子间的真实剪切速度差施加粘性阻尼，排除刚体旋转带来的切向速度差，同时保留径向速度差的耗散（满足行星聚集需求）。以下是两种高效实现方案：

方案1：单次遍历的剪切分量过滤法 #

无需额外循环，在现有邻居遍历流程中即可完成计算。核心是分离相对速度中的刚体旋转分量（无剪切变形）与真实剪切分量：

// 遍历单个邻居粒子时的计算逻辑
vec3 r_ij = pos_neighbor - pos_particle;
float r = length(r_ij);
if (r < 1e-6) continue; // 避免除以0
vec3 r_hat = r_ij / r;
vec3 v_ij = vel_neighbor - vel_particle;

// 分离相对速度的剪切分量：移除径向分量后剩余的部分
vec3 v_shear = v_ij - dot(v_ij, r_hat) * r_hat;

// 仅对剪切分量施加粘性力
visc_force += viscosity * mass_neighbor / density_neighbor * v_shear * kernelFun_laplacianf(r, r0);

效果说明

• 刚体旋转时，相邻粒子的相对速度完全由刚体运动产生，无真实剪切变形，v_shear为0，粘性力不生效，不会阻尼旋转；
• 层流等真实剪切流动中，v_shear对应实际剪切速度差，粘性力会正确耗散剪切能量，维持流体行为；
• 行星聚集时的径向压缩/膨胀速度差被排除在v_shear外，你可以选择保留径向分量的粘性耗散（只需在代码中叠加径向分量的贡献），完全满足行星形成的能量需求。

方案2：基于应变率张量的精确粘性计算 #

如果需要更精确的流体粘性行为（如符合Navier-Stokes方程的粘性），可以基于应变率张量计算，同样可在单次邻居遍历中完成（需累积速度梯度）：

// 第一步：遍历邻居累积速度梯度
mat3 vel_grad = mat3(0.0);
for (每个邻居j) {
    vec3 r_ij = pos_neighbor - pos_particle;
    float r = length(r_ij);
    vec3 grad_W = kernelFun_gradient(r, r0); // 核函数的梯度
    vec3 v_ij = vel_neighbor - vel_particle;
    vel_grad += (mass_neighbor / density_neighbor) * outer_product(v_ij, grad_W);
}

// 计算应变率张量（对称速度梯度，仅包含剪切和膨胀变形）
mat3 strain_rate = 0.5 * (vel_grad + transpose(vel_grad));

// 第二步：遍历邻居计算粘性力（可与第一步合并为单次遍历，预存邻居信息）
vec3 visc_force = vec3(0.0);
for (每个邻居j) {
    vec3 r_ij = pos_neighbor - pos_particle;
    vec3 grad_W = kernelFun_gradient(r, r0);
    float rho_avg = 0.5 * (density_particle + density_neighbor);
    // 粘性应力张量的贡献
    vec3 stress_contrib = 2.0 * viscosity * strain_rate * r_ij;
    visc_force += (mass_neighbor / rho_avg) * stress_contrib * grad_W;
}

效果说明

应变率张量仅包含流体的变形部分，自动排除刚体旋转的反对称速度梯度分量，因此粘性力只会作用于真实的剪切和膨胀变形，完全符合流体力学规律，同时避免刚体旋转阻尼。

方案对比 #

• 方案1实现简单，计算量极小，完全适配现有遍历逻辑，无额外开销；
• 方案2更精确，符合Navier-Stokes方程，适合对流体行为要求较高的场景；
• 两者均不会破坏层流等真实流体效果，同时解决刚体旋转阻尼问题，满足行星聚集的能量耗散需求。

内容的提问来源于stack exchange，提问作者Paul Aner