SPH（Smoothed Particle Hydrodynamics，平滑粒子流体动力学）是一种无网格的拉格朗日方法，广泛用于流体动力学分析。
1. 理解SPH方法的基本原理
- 粒子表示流体：流体被离散为一系列粒子，每个粒子携带质量、位置、速度、密度等属性。
- 平滑函数：通过平滑核函数（如高斯核函数）对粒子的物理量进行插值，从而计算流体的物理性质。
- 运动方程：基于牛顿第二定律和流体动力学方程，更新粒子的位置和速度。
2. MATLAB实现SPH方法的步骤
（1）初始化粒子
- 定义流体粒子的初始位置、速度、密度等属性。
% 参数设置
numParticles = 1000; % 粒子数量
domainSize = [0, 1; 0, 1]; % 流体域范围
particleMass = 1e-4; % 每个粒子的质量
h = 0.05; % 平滑长度
rho0 = 1000; % 参考密度（水的密度）
gamma = 7; % 状态方程参数
cs = 10; % 声速
% 初始化粒子位置和速度
positions = zeros(numParticles, 2);
velocities = zeros(numParticles, 2);
% 随机生成粒子位置
for i = 1:numParticles
positions(i, :) = domainSize(:, 1) + (domainSize(:, 2) - domainSize(:, 1)) .* rand(2, 1);
end
（2）定义平滑核函数
- 平滑核函数用于计算粒子之间的相互作用。
% 高斯平滑核函数
function W = gaussianKernel(r, h)
sigma = h / 1.2;
W = (1 / (2 * pi * sigma^2)) * exp(-r^2 / (2 * sigma^2));
end
% 核函数的梯度
function gradW = gradGaussianKernel(r, h)
sigma = h / 1.2;
gradW = (-r / sigma^2) * gaussianKernel(r, h);
end
（3）计算粒子密度和压力
- 使用平滑核函数计算每个粒子的密度。
- 根据状态方程计算压力。
% 计算密度
densities = zeros(numParticles, 1);
for i = 1:numParticles
for j = 1:numParticles
if i ~= j
r = norm(positions(i, :) - positions(j, :));
densities(i) = densities(i) + particleMass * gaussianKernel(r, h);
end
end
end
% 计算压力
pressures = zeros(numParticles, 1);
for i = 1:numParticles
pressures(i) = cs^2 * ((densities(i) / rho0)^(gamma) - 1);
end
（4）计算粒子间相互作用力
- 根据压力梯度计算粒子间的相互作用力。
% 计算相互作用力
forces = zeros(numParticles, 2);
for i = 1:numParticles
for j = 1:numParticles
if i ~= j
r = positions(i, :) - positions(j, :);
dist = norm(r);
gradW = gradGaussianKernel(dist, h);
pressureForce = -particleMass * (pressures(i) / densities(i)^2 + pressures(j) / densities(j)^2) * gradW * r / dist;
forces(i, :) = forces(i, :) + pressureForce;
end
end
end
（5）更新粒子位置和速度
- 使用欧拉方法或其他时间积分方法更新粒子的位置和速度。
% 时间步长
dt = 1e-4;
% 更新速度
velocities = velocities + forces / densities * dt;
% 更新位置
positions = positions + velocities * dt;
（6）可视化结果
- 使用 MATLAB 的绘图功能显示流体粒子的分布。
figure;
scatter(positions(:, 1), positions(:, 2), 36, 'filled');
title('SPH 流体模拟');
xlabel('X');
ylabel('Y');
axis equal;
3. 完整代码示例
以下是将上述步骤整合到一个完整的 MATLAB 脚本中的示例代码：
% 参数设置
numParticles = 1000; % 粒子数量
domainSize = [0, 1; 0, 1]; % 流体域范围
particleMass = 1e-4; % 每个粒子的质量
h = 0.05; % 平滑长度
rho0 = 1000; % 参考密度（水的密度）
gamma = 7; % 状态方程参数
cs = 10; % 声速
dt = 1e-4; % 时间步长
numSteps = 100; % 模拟步数
% 初始化粒子位置和速度
positions = zeros(numParticles, 2);
velocities = zeros(numParticles, 2);
for i = 1:numParticles
positions(i, :) = domainSize(:, 1) + (domainSize(:, 2) - domainSize(:, 1)) .* rand(2, 1);
end
% 平滑核函数
function W = gaussianKernel(r, h)
sigma = h / 1.2;
W = (1 / (2 * pi * sigma^2)) * exp(-r^2 / (2 * sigma^2));
end
% 核函数的梯度
function gradW = gradGaussianKernel(r, h)
sigma = h / 1.2;
gradW = (-r / sigma^2) * gaussianKernel(r, h);
end
% 主循环
for step = 1:numSteps
% 计算密度
densities = zeros(numParticles, 1);
for i = 1:numParticles
for j = 1:numParticles
if i ~= j
r = norm(positions(i, :) - positions(j, :));
densities(i) = densities(i) + particleMass * gaussianKernel(r, h);
end
end
end
% 计算压力
pressures = zeros(numParticles, 1);
for i = 1:numParticles
pressures(i) = cs^2 * ((densities(i) / rho0)^(gamma) - 1);
end
% 计算相互作用力
forces = zeros(numParticles, 2);
for i = 1:numParticles
for j = 1:numParticles
if i ~= j
r = positions(i, :) - positions(j, :);
dist = norm(r);
gradW = gradGaussianKernel(dist, h);
pressureForce = -particleMass * (pressures(i) / densities(i)^2 + pressures(j) / densities(j)^2) * gradW * r / dist;
forces(i, :) = forces(i, :) + pressureForce;
end
end
end
% 更新速度和位置
velocities = velocities + forces / densities * dt;
positions = positions + velocities * dt;
% 可视化
if mod(step, 10) == 0
figure;
scatter(positions(:, 1), positions(:, 2), 36, 'filled');
title(['SPH 流体模拟 - 步骤 ', num2str(step)]);
xlabel('X');
ylabel('Y');
axis equal;
end
end
参考代码 采用matlab，SPH方法进行流体分析 www.youwenfan.com/contentcno/53224.html
4. 注意事项
- 平滑长度 (h) 的选择：平滑长度 (h) 影响粒子之间的相互作用范围，需要根据问题规模合理选择。
- 时间步长 (dt) 的选择：时间步长需要足够小以保证数值稳定性，通常需要满足 CFL 条件。
- 边界处理：在实际应用中，需要考虑流体与边界的相互作用，可以通过添加虚拟粒子或施加边界条件来实现。
