哪位大侠有sph的源程序共享一下

mechanics

刚开始接触，哪位大侠有源程序吗，可否共享一下

sdfluid

简单地说，SPH方法就是用粒子离散计算域的场变量。那么，所有场变量（密度、速度、压力等）就都必须附到每一个SPH粒子上，成为该粒子的属性值。而该属性值是以包含该粒子的一个局部区域内的其他粒子的对应属性值加权求和得到的。
Cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
C                               SPH算法 第一部分
Cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
所以，SPH算法的第一步骤是局部区域（影响域）内粒子搜索，书中称之“nearest neighboring particle search（NNPS）”。
对计算域内所有粒子i都进行NNPS搜索，可采取成对相互作用法降低工作量。
do i=1,ntotal-1
  do j=i+1,ntotal
通过距离比较判定任意一个粒子j是否在粒子i的影响域内。
If ( |x(i)-x(j)|)<=k_factor*h)
Niac++;
//  niac总数=（ntotal-1,…ntotal-(ntotal-1)降序排列求和=n(n-1)/2）
// 粒子数过千，碰撞对数上十万。而且每个时间步都还要用此碰撞对计算一遍。
因为后文需SPH近似N-S方程，故这里要把所有碰撞对存储下。
   Pair_i(niac)=I;
   Pair_j(niac)=j;
//matlab试验下，这个设计非常妙。
// count=0;
n=10;
p(1:45)=0;
q(1:45)=0;
for i=1:n-1;
    for j=i+1:n;
        count=count+1;
        p(count)=i;
        q(count)=j;
    end
end
// p= [1     1     1     1     1     1     1     1     1     2     2     2     2     2     2     2     2     3     3     3     3     3     3     3     4   4     4     4     4     4     5     5     5     5     5     6     6     6     6     7     7     7     8     8     9]
// q= [2     3     4     5     6     7     8     9    10     3     4     5     6     7     8     9    10     4     5     6     7     8     9    10     5   6     7     8     9    10     6     7     8     9    10     7     8     9    10     8     9    10     9    10    10 ]
// SPH算法第二部分也大量调用了这个结构。
每个碰撞对都调用一次核函数，计算Wij和导数，为后文SPH近似做准备
   Call kernel(..)
   Endif
  Enddo
//内层do完成，得到将对一个特定i粒子属性值贡献的所有J粒子
Enddo
//外层do完成，相当于确认了全域每个粒子的影响域内的其他粒子。各粒子在空间的位置还没变，相当于全域在时间轴上位于（当前时间步的）初始时刻t0

C******************************************************************************
C                           SPH 算法 第二部分
C******************************************************************************
SPH算法的第一步，完成了将场变量转换成粒子属性。第二步，是一个时间步内的场变量更新。
变量更新的控制方程就是N-S方程中的动量守恒方程。该方程涉及到的场变量包括：密度、压力、速度（用于计算剪应变）以及粘性系数和壁面斥力。
书中将N-S方程右端项，即压力梯度、物理粘性力称为内力（Internal force）；壁面斥力（以及重力场）称为外力(external force)。N-S方程左端即加速度项。

首先，SPH近似密度场
Do i=1,ntotal
Call kernel(r,dr,h,selfdens,hdwdx)
Rho(i)=selfdens*mass(i)
Enddo
// 直接令rho(i)=0，会有多大精度损失？
Do k=1,niac
I=pair_i(k)
J=pair_j(k)
Rho(i)=rho(i)+mass(j)*w(k)
Rho(j)=rho(j)+mass(i)*w(k)
enddo
//MATLAB验算
// count=0;
n=10;
p(1:45)=0;
q(1:45)=0;
for i=1:n-1;
    for j=i+1:n;
        count=count+1;
        p(count)=i;
        q(count)=j;
    end
end
r(1:45)=0;
for k=1:45;
    i=p(k);
    j=q(k);
    r(i)=r(i)+j;
    r(j)=r(j)+i;
    r(i)
    r(j)
end
R(i)=[2,5,9,14,20,27,35,44,54;4,8,13,19,26,43,53;7,12,18,25,33,42,52;11,17,24,32,41,51;
16,23,31,40,50;22,30,39,49;29,38,48;37,47;46]
含义表述： do k=1:9, r(i)前九个数对应与所有与i粒子发生碰撞的其他粒子对i粒子的属性r（i）贡献累加。
以密度的SPH近似为例，其实现了:i处粒子密度有影响域内其他粒子加权求和表示。
同理其他场变量。后文不重复讨论。

//下面，计算NS方程右端项贡献的加速度
do k=1,niac
I=pair_i(k)
J=pair_j(k)
c.. calculate strain for viscosity force
c.. calculate pressure force
return indvdxt(d,i)
enddo
//下面计算外力项对加速度的贡献
Do i=1,ntotal
Dvxdt(dim,i)=-9.8
enddo
do k=1,niac
I=pair_i(k)
J=pair_j(k)
c.. calculate boundary repulse force
return exdvdxt(d,i)
enddo
//加速度贡献各项求和，并进行时间积分，求速度项
Do i=1,ntotal
   Do d=1,dim
   Dvx(d,i)=indvdxt(d,i)+exdvdxt(d,i)
   Enddo
Enddo
// 时间积分，求当前时间步末尾时刻(t+dt)各粒子的速度、位移
  Do i=1,ntotal
     Do d=1,dim
     Vx(d,i)=v0(d,i)+dt*dvx(d,i)
     X(d,i)=x0(d,i)+dt*vx(d,i)
    Enddo
  Enddo
//进入下一个时间步
  Time=time+dt

sara

sdfluid 发表于 2013-11-15 11:09
简单地说，SPH方法就是用粒子离散计算域的场变量。那么，所有场变量（密度、速度、压力等）就都必须附到每一 ...

算什么的这个

wonex

2楼不错，研究一下