如何让SPH粒子更符合流体特性？JavaScript流体模拟优化问询

你的模拟出现粒子随机排斥、加重力后混乱的问题，核心是公式实现错误、参数配置不合理和初始化逻辑问题，以下是针对性修复：

1. 密度计算的self引用错误 #

calculatedensity方法中用self判断粒子是否为自身，但JavaScript类方法里必须用this，导致粒子计算密度时把自己也纳入计算，密度值异常。

修复代码：

calculatedensity(nearbyparticles) {
  let density = 0;
  for (let cparticle of nearbyparticles) {
    // 用this代替self
    if (cparticle !== this) {
      let distance = this.position.distanceFrom(cparticle.position);
      density += cparticle.mass * fluidSimulator.poly6kernel(distance, this.smoothinglength);
    }
  }
  this.density = density;
}

2. 压力力公式错误 #

SPH的压力力需要同时考虑当前粒子和邻居粒子的压力（压力是相互作用），原代码仅使用邻居粒子的压力，且核函数的应用逻辑有误，导致力的方向和大小异常。

修复代码：

calculatepressureforce(nearbyparticles) {
  let pressureforce = new Vector(0, 0);
  for (let cparticle of nearbyparticles) {
    if (cparticle === this) continue;
    let distance = this.position.distanceFrom(cparticle.position);
    if (distance > this.smoothinglength || distance === 0) continue;

    // 计算从邻居指向当前粒子的方向向量
    let direction = cparticle.position.sub(this.position).divScalar(distance);
    // 正确的压力力公式：平均压力 * 质量 * 核导数 / 邻居密度
    let avgPressure = (this.pressure + cparticle.pressure) / 2;
    let magnitude = avgPressure * cparticle.mass * fluidSimulator.spikykernelderrivative(distance, this.smoothinglength) / cparticle.density;
    
    pressureforce = pressureforce.add(direction.mulScalar(magnitude));
  }
  // 压力力是排斥力，无需额外取负（核导数已包含符号）
  this.pressureforce = pressureforce;
}

3. 参数配置优化 #

原参数导致力的强度和影响范围不合理，调整后保证数值稳定性：

• 光滑长度改为3*radius（SPH标准取值，避免过度影响）
• 压力乘数改为2（避免压力过强）
• 碰撞阻尼改为0.5（碰撞后损失部分能量，避免弹开过度）
• 时间步长改为1/60（更小的步长提升稳定性）
• 修正逻辑与运算符：把&改为&&，避免位运算错误

修改参数部分：

var radius = 0.2;
var pressuremultiplier = 2; // 从10改为2
var restdensity = 1;

// Particle类构造函数中
this.smoothinglength = 3 * radius; // 从6*radius改为3*radius

// SPHFluidSimulator构造函数中
this.timestep = 1 / 60; // 从1/30改为1/60
this.collisionDamping = 0.5; // 从1改为0.5

4. 粒子初始化改为网格分布 #

随机初始化的粒子初始密度不均匀，导致初始压力爆发，改用网格排列让粒子初始均匀分布：

替换generateParticlesRandom函数：

function generateParticlesGrid() {
  let gridSize = Math.sqrt(numParticles);
  let spacing = 0.3; // 略大于粒子半径，保证初始均匀
  for (let i = 0; i < numParticles; i++) {
    let x = (i % gridSize) * spacing + spacing;
    let y = Math.floor(i / gridSize) * spacing + spacing;
    // 限制在模拟区域内
    x = Math.min(x, simWidth - spacing);
    y = Math.min(y, simHeight - spacing);
    let particle = new Particle(x, y);
    fluidSimulator.initialiseParticle(particle);
  }
}

// 调用时替换generateParticlesRandom()
generateParticlesGrid();

完整修复代码 #

var canvas = document.getElementById("canvas");
var c = canvas.getContext("2d");
canvas.width = window.innerWidth - 20;
canvas.height = window.innerHeight - 100;
var simMinWidth = 20.0;
var cScale = Math.min(canvas.width, canvas.height) / simMinWidth;
var simWidth = canvas.width / cScale;
var simHeight = canvas.height / cScale;
var numParticles = 100;
var radius = 0.2;
var pressuremultiplier = 2;
var restdensity = 1;

function cX(position) {
  return position.x * cScale;
}

function cY(position) {
  return position.y * cScale;
}

class Particle {
  constructor(x, y) {
    this.position = new Vector(x, y);
    this.velocity = new Vector(0, 0);
    this.radius = radius;
    this.smoothinglength = 3 * radius;
    this.density = 1;
    this.pressure = 1;
    this.pressureforce = new Vector(0, 0); // 初始化为向量而非数值
    this.mass = 1;
  }
  updateAcceleration() {
    this.calculatedensity(fluidSimulator.particles);
    this.calculatepressure();
    this.calculatepressureforce(fluidSimulator.particles);
    this.acceleration = this.pressureforce.divScalar(this.density).add(fluidSimulator.gravity);
  }
  checkboundarys() {
    if (this.position.x - this.radius < 0.0) {
      this.position.x = 0.0 + this.radius;
      this.velocity.x = -fluidSimulator.collisionDamping * this.velocity.x;
    }

    if (this.position.x + this.radius > simWidth) {
      this.position.x = simWidth - this.radius;
      this.velocity.x = -fluidSimulator.collisionDamping * this.velocity.x;
    }

    if (this.position.y - this.radius < 0.0) {
      this.position.y = 0.0 + this.radius;
      this.velocity.y = -fluidSimulator.collisionDamping * this.velocity.y;
    }

    if (this.position.y + this.radius > simHeight) {
      this.position.y = simHeight - this.radius;
      this.velocity.y = -fluidSimulator.collisionDamping * this.velocity.y;
    }
  }
  calculatedensity(nearbyparticles) {
    let density = 0;
    for (let cparticle of nearbyparticles) {
      if (cparticle !== this) {
        let distance = this.position.distanceFrom(cparticle.position);
        density += cparticle.mass * fluidSimulator.poly6kernel(distance, this.smoothinglength);
      }
    }
    this.density = density;
  }
  calculatepressure() {
    this.pressure = pressuremultiplier * (this.density - restdensity);
  }
  calculatepressureforce(nearbyparticles) {
    let pressureforce = new Vector(0, 0);
    for (let cparticle of nearbyparticles) {
      if (cparticle === this) continue;
      let distance = this.position.distanceFrom(cparticle.position);
      if (distance > this.smoothinglength || distance === 0) continue;

      let direction = cparticle.position.sub(this.position).divScalar(distance);
      let avgPressure = (this.pressure + cparticle.pressure) / 2;
      let magnitude = avgPressure * cparticle.mass * fluidSimulator.spikykernelderrivative(distance, this.smoothinglength) / cparticle.density;
      
      pressureforce = pressureforce.add(direction.mulScalar(magnitude));
    }
    this.pressureforce = pressureforce;
  }

  simulateParticle() {
    this.updateAcceleration();
    this.velocity = this.velocity.add(this.acceleration.mulScalar(fluidSimulator.timestep));
    this.position = this.position.add(this.velocity.mulScalar(fluidSimulator.timestep));
    this.checkboundarys();
  }
}

class SPHFluidSimulator {
  constructor() {
    this.timestep = 1 / 60;
    this.gravity = new Vector(0, 0); // 可改为new Vector(0, 9.8)添加重力
    this.collisionDamping = 0.5;
    this.particleSpacing = 0.1;
    this.particles = [];
  }
  poly6kernel(r, h) {
    let influence = 0;
    if (r >= 0 && r <= h) { // 用&&代替&，避免位运算错误
      influence = Math.pow(h ** 2 - r ** 2, 3);
    }
    return (315 / (64 * Math.PI * Math.pow(h, 9))) * influence;
  }
  spikykernelderrivative(r, h) {
    let influence = 0;
    if (r >= 0 && r <= h) { // 用&&代替&
      influence = -3 * Math.pow(h - r, 2);
    }
    return (15 / (Math.PI * Math.pow(h, 6))) * influence;
  }
  initialiseParticle(particle) {
    this.particles.push(particle);
  }
}

function generateParticlesGrid() {
  let gridSize = Math.sqrt(numParticles);
  let spacing = 0.3;
  for (let i = 0; i < numParticles; i++) {
    let x = (i % gridSize) * spacing + spacing;
    let y = Math.floor(i / gridSize) * spacing + spacing;
    x = Math.min(x, simWidth - spacing);
    y = Math.min(y, simHeight - spacing);
    let particle = new Particle(x, y);
    fluidSimulator.initialiseParticle(particle);
  }
}

function physics() {
  for (const particle of fluidSimulator.particles) {
    particle.simulateParticle();
  }
}

function draw() {
  c.clearRect(0, 0, canvas.width, canvas.height);
  c.fillStyle = "#005EB8";
  for (let i = 0; i < numParticles; i++) {
    let particle = fluidSimulator.particles[i];
    c.beginPath();
    c.arc(cX(particle.position), cY(particle.position), cScale * radius, 0.0, 2.0 * Math.PI);
    c.closePath();
    c.fill();
  };
}

function simulate() {
  physics();
  draw();
  requestAnimationFrame(simulate);
}

var fluidSimulator = new SPHFluidSimulator();
generateParticlesGrid();
// 添加重力可取消注释：
// fluidSimulator.gravity = new Vector(0, 9.8);
simulate();

<script src="https://ronenness.github.io/Vector2js/vector2js.js"></script>
<canvas id="canvas"></canvas>

效果说明 #

• 无重力时：粒子会逐渐调整到均匀分布状态，无聚集或大间隙
• 添加重力（设置fluidSimulator.gravity = new Vector(0,9.8)）：粒子会自然下落，碰撞边界后阻尼反弹，呈现流体流动的效果，不会出现爆炸式混乱

内容的提问来源于stack exchange，提问作者Splashking XD