1
2
3
4
5
6
7
8
9

struct SPHParticle
{
    float3 position;    // 位置
    float3 velocity;    // 速度
    float3 force;       // 受力
    float density;      // 密度
    float pressure;     // 压力
    float2 grid;        // 网格坐标（用于后续优化）
};

1
2
3

1. 初始化粒子位置 → 2. 计算密度和压力 → 3. 计算受力
    ↓                                    ↓
4. 更新速度和位置 ← 5. 碰撞检测和响应

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

 Assets/
 ├── Script/
 │   └── MySph/
 │       ├── Base/                    # 基础数据结构
 │       │   └── SphParticle.cs
 │       ├── SPHSolver/
 │       │   └── SphSolver.cs          # 求解器主类
 │       ├── SphParticleSystem.cs      # 粒子系统管理
 │       └── SphColliderMesh.cs        # 碰撞体组件
 ├── Resources/
 │   └── ComputeShader/
 │       └── Sph/
 │           └── SphSolver.compute     # 计算着色器
 └── Meshes/                           
     └── SphereLowPoly.fbx             # 粒子 Mesh

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14

// 获取粒子系统组件
var particleSystem = GetComponent<SphParticleSystem>();

// 配置粒子数量和半径
particleSystem.particleCount = 8000;
particleSystem.smoothingRadius = 1.0f;
particleSystem.particleMass = 0.1f;

// 添加碰撞体
var collider = gameObject.AddComponent<SphColliderMesh>();
collider.transform.Find("CollidableMesh").GetComponent<MeshCollider>();

// 添加求解器
var solver = gameObject.AddComponent<SphSolver>();

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74

float Wpoly6(float3 r, float h)
{
    return 315.0f / (64.0f * PI * pow(h, 9.0f)) * pow(h * h - dot(r, r), 3.0f) * step(length(r), h);
}

[numthreads(PARTICLETHREADCOUNT,1,1)]
void ComputeDensityPressure(uint3 id : SV_DispatchThreadID)
{
    SPHParticle particle = particles[id.x];
    particle.density = 0.0f;

    for (int j = 0; j < particleCount; j++)
    {
        particle.density += mass * Wpoly6(particles[j].position - particle.position, smoothingRadius);
    }
    // 计算pressure
    particle.pressure = gas * (particle.density - restDensity);

    particles[id.x] = particle;
}


float3 SinglePressureForce(float3 ri, float3 rj, float pi, float pj, float di, float dj)
{
    float r = length(ri - rj);
    return (pi + pj) / (2.0f * dj) * pow(smoothingRadius - r, 2.0f) * (ri - rj) / r * step(r, smoothingRadius);
}

float3 ComputePressureForce(uint i)
{
    float3 forcePressure = 0;
    for (uint j = 0; j < particleCount; j++)
        forcePressure += i == j
                             ? 0
                             : SinglePressureForce(
                                 particles[i].position, particles[j].position,
                                 particles[i].pressure, particles[j].pressure,
                                 particles[i].density, particles[j].density);
    // Spiky 核函数
    forcePressure *= mass * (-45.0f / (PI * pow(smoothingRadius, 6.0f)));
    return forcePressure;
}

float3 SingleViscosityForce(float3 ri, float3 rj, float di, float dj, float vi, float vj)
{
    float r = length(ri - rj);
    return (smoothingRadius - r) * (vj - vi) / dj * step(r, smoothingRadius);
}

float3 ComputeViscosityForce(uint i)
{
    float3 forceViscosity = 0;
    for (uint j = 0; j < particleCount; j++)
        forceViscosity += i == j
                              ? 0
                              : SingleViscosityForce(
                                  particles[i].position,
                                  particles[j].position,
                                  particles[i].density,
                                  particles[j].density,
                                  particles[i].velocity,
                                  particles[j].velocity);
    // ∇²W_viscosity 核函数
    forceViscosity *= mass * particleViscosity * (45.0f / (PI * pow(smoothingRadius, 6.0f)));
    return forceViscosity;
}

[numthreads(PARTICLETHREADCOUNT,1,1)]
void ComputeForces(uint3 id : SV_DispatchThreadID)
{
    float3 forceGravity = gravity.xyz * particles[id.x].density * gravity.w;
    particles[id.x].force = ComputePressureForce(id.x) + ComputeViscosityForce(id.x) + forceGravity;
}