📄 sph.cpp
字号:
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
#include "SPH.h"
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
SPH::SPH(int dim1,int maxn1)
{
dim=dim1;
maxn=maxn1;
max_interaction = 100 * maxn;
x_maxgeom =10.0;
x_mingeom = -10.0;
y_maxgeom =10.0;
y_mingeom = -10.0;
z_maxgeom =10.0;
z_mingeom = -10.0;
pa_sph = 2;
nnps = 1;
sle = 0;
skf = 1;
summation_density = true;
average_velocity = false;
config_input = false;
virtual_part = false;
vp_input = false;
visc = false;
ex_force =false;
visc_artificial = true;
heat_artificial = false;
self_gravity = false;
nor_density = false;
nsym = 0;
int_stat = true;
print_step = 100;
save_step = 500;
moni_particle = 1600;
pi = 3.14159265358979323846;
shocktube = true;
shearcavity = false;
itimestep=0;//新添加
//----------------input--------------------
p=new double[maxn+1];
u=new double[maxn+1];
//----------------direct_find--------------------
countiac=new int [maxn+1];
//--------------kernel-----------------
dx=new double[dim+1];
// dw1dx=new double[dim+1];
//------art_heat----------------------
// ntotal= ;//?????????????有子程序定义，此处不必再定义
// niac= ;//??????????????有子程序定义，此处不必再定义
pair_i=new int[max_interaction+1];
pair_j=new int[max_interaction+1];
hsml=new double[maxn+1];
mass=new double[maxn+1];
x=new double *[dim+1];
for(int i=0;i<dim+1;i++)
x[i]=new double[maxn+1];
vx=new double *[dim+1];
for(i=0;i<dim+1;i++)
vx[i]=new double[maxn+1];
rho=new double[maxn+1];
u=new double[maxn+1];
c=new double[maxn+1];
w=new double[max_interaction+1];
dwdx=new double *[dim+1];
for(i=0;i<dim+1;i++)
dwdx[i]=new double[max_interaction+1];
// dedt=new double[maxn];
//---------------art_vsic------------------
// dvxdt=new double *[dim+1];
// for(i=0;i<dim+1;i++)
// dvxdt[i]=new double[maxn+1];
//----------------av_vel---------------------
av=new double *[dim+1];
for(i=0;i<dim+1;i++)
av[i]=new double[maxn+1];
//---------------------sum_density------------
itype=new int[maxn+1];
//---------------------con_density------------
drhodt=new double[maxn+1];
//-----------------------------------int_force--------------------------------------
//此程序中出现两个itype的定义，请查明原因
eta=new double[maxn+1];
t=new double[maxn+1];
tdsdt=new double[maxn+1];
//-----------------------------------single_step--------------------------------------
s=new double[maxn+1];//此变量本程序中没有使用
du=new double[maxn+1];
ds=new double[maxn+1];
drho=new double[maxn+1];
dvx=new double *[dim+1];
for(i=0;i<dim+1;i++)
dvx[i]=new double[maxn+1];
//-----------------------------------time_integration--------------------------------------
//e=new double[maxn+1];
}
SPH::~SPH()
{
//-----------------------------------int_force--------------------------------------
delete []eta;
delete []t;
delete []tdsdt;
//------------------------------------------------------------------------------
//-----------------------------------single_step--------------------------------------
delete []s;//此变量本程序中没有使用
delete []du;
delete []ds;
delete []drho;
for(int i=0;i<dim+1;i++)
delete []dvx[i];
delete []dvx;
//------------------------------------------------------------------------------
//-----------------------------------time_integration--------------------------------------
//delete []e;
//--------------input-----------
delete []p;
// delete []u; //重复定义
//--------------direct_find-----------
delete []countiac;
//---------------kernel---------------
delete[]dx;
// delete[]dw1dx;
//------------------------------------
//------art_heat----------------------
delete []pair_i;
delete []pair_j;
delete []hsml;
delete []mass;
for( i=0;i<dim+1;i++)
delete []x[i];
delete []x;
for( i=0;i<dim+1;i++)
delete []vx[i];
delete []vx;
delete []rho;
delete []u;
delete []c;
delete []w;
for( i=0;i<dim+1;i++)
delete []dwdx[i];
delete []dwdx;
// delete []dedt;
//-------------------------------------
//------------av_vel-------------------
for( i=0;i<dim+1;i++)
delete []av[i];
delete []av;
//--------------------------------------
//--------------sum_density-------------
delete []itype;
//------------------------------------
//--------------con_density-------------
delete []drhodt;
//------------------------------------
//---------------art_vsic------------------
// for(i=0;i<dim+1;i++)
// delete []dvxdt[i];
// delete []dvxdt;
//-------------------------------------------
}
void p_gas(double rho1,double u1,double &p1,double &c1)
{
double gamma;
gamma=1.4;
p1 = (gamma-1) * rho1 * u1;
c1 = sqrt((gamma-1) * u1);
}
void p_art_water(double rho1,double &p1,double &c1)
{
//double gamma=7;
//double rho0=1000;
c1 = 0.01;
p1 = pow(c1,2) * rho1;
}
void SPH::art_visc(int ntotal,double **(&dvx1dt),double *(&de1dt))
{
// if (visc_artificial)art_visc(ntotal+nvirt,hsml,mass,x,vx,niac,rho,c,pair_i,pair_j,w,dwdx,ardvxdt,avdudt);//其中c表示温度还是声速，这里以声速计算
//c Subroutine to calculate the artificial viscosity (Monaghan, 1992)
// ntotal : Number of particles (including virtual particles) [in]
// hsml : Smoothing Length [in]
// mass : Particle masses [in]
// x : Coordinates of all particles [in]
// vx : Velocities of all particles [in]
// niac : Number of interaction pairs [in]
// rho : Density [in]
// c : Temperature [in]
// pair_i : List of first partner of interaction pair [in]
// pair_j : List of second partner of interaction pair [in]
// w : Kernel for all interaction pairs [in]
// dwdx : Derivative of kernel with respect to x, y and z [in]
// dvxdt : Acceleration with respect to x, y and z [out]
// dedt : Change of specific internal energy [out]
double dx1;//已定义*dx,故改名
double piv,muv, vr, rr, h, mc, mrho, mhsml;
double *dvx1;// 构造函数dvx已定为二维数组，故改名，dvx1[dim]不能直接定义数组是因为数组的大小在构造函数并没有指明大小（也是动态的）
dvx1=new double[dim+1];
// Parameter for the artificial viscosity:
double alpha = 1.0;//Shear viscosity查明这三个参数是不是在其他程序也引用，这是当作局部变量
double beta = 1.0; // Bulk viscosity
double etq = 0.1;// Parameter to avoid singularities
for(int i=1;i<=ntotal;i++)
{
for(int d=1;d<=dim;d++)
{
dvx1dt[d][i] = 0.0;
}
de1dt[i] = 0.0;
}
// Calculate SPH sum for artificial viscosity
for(int k=1;k<=niac;k++)
{
i = pair_i[k];
int j = pair_j[k];
mhsml= (hsml[i]+hsml[j])/2;
vr = 0.0;
rr = 0.0;
for(int d=1;d<=dim;d++)
{
dvx1[d]=vx[d][i]-vx[d][j];
dx1=x[d][i]-x[d][j];
vr=vr+dvx1[d]*dx1;
rr=rr+dx1*dx1;
}
// Artificial viscous force only if v_ij * r_ij < 0
if (vr<0.0)
{
// Calculate muv_ij = hsml v_ij * r_ij / ( r_ij^2 + hsml^2 etq^2 )
muv = mhsml*vr/(rr + mhsml*mhsml*etq*etq);
// Calculate PIv_ij = (-alpha muv_ij c_ij + beta muv_ij^2) / rho_ij
mc = 0.5*(c[i] + c[j]);
mrho = 0.5*(rho[i] + rho[j]);
piv = (beta*muv - alpha*mc)*muv/mrho;
// Calculate SPH sum for artificial viscous force
for(int d=1;d<=dim;d++)
{
h = -piv*dwdx[d][k];
dvx1dt[d][i] = dvx1dt[d][i] + mass[j]*h;
dvx1dt[d][j] = dvx1dt[d][j] - mass[i]*h;
de1dt[i] = de1dt[i] - mass[j]*dvx1[d]*h;
de1dt[j] = de1dt[j] - mass[i]*dvx1[d]*h;
}
}
}
// Change of specific internal energy:
for(i=1;i<=ntotal;i++)
de1dt[i] = 0.5*de1dt[i];
// for(i=1;i<=ntotal;i++)cout<<"de1dt[i]="<<de1dt[i]<<"dvx1dt[1][i]="<<dvx1dt[1][i]<<endl;
delete []dvx1;
}
void SPH::sum_density(int ntotal)
{
// sum_density(ntotal+nvirt,hsml,mass,niac,pair_i,pair_j,w,itype,rho);
// C Subroutine to calculate the density with SPH summation algorithm.
//C ntotal : Number of particles [in]
//C hsml : Smoothing Length [in]
//C mass : Particle masses [in]
//C niac : Number of interaction pairs [in]
//C pair_i : List of first partner of interaction pair [in]
//C pair_j : List of second partner of interaction pair [in]
//C w : Kernel for all interaction pairs [in]
//c itype : type of particles [in]
//c x : Coordinates of all particles [in]
//c rho : Density [out]
double selfdens,r;
double *hv, *wi;
hv=new double[dim+1];
wi=new double[maxn+1];
// wi[maxn]---integration of the kernel itself
for(int d=1;d<=dim;d++)
hv[d] = 0.0;//这个变量程序中好像并没有使用？？？？？？？？？？
// Self density of each particle: Wii (Kernel for distance 0)
// and take contribution of particle itself:
r=0;
// Firstly calculate the integration of the kernel over the space
double temp=0;
for(int i=1;i<=ntotal;i++)
{
temp=hsml[i];
kernel(r,hv,temp,selfdens,hv);//这是程序可能有问题的一个地方??????????????????????????????
wi[i]=selfdens*mass[i]/rho[i];
}
for(int k=1;k<=niac;k++)
{
i = pair_i[k];
int j = pair_j[k];
wi[i] = wi[i] + mass[j]/rho[j]*w[k];
wi[j] = wi[j] + mass[i]/rho[i]*w[k];
}
// Secondly calculate the rho integration over the space
for( i=1;i<=ntotal;i++)
{
temp=hsml[i];
kernel(r,hv,temp,selfdens,hv);//这是程序可能有问题的一个地方???????????????????
rho[i] = selfdens*mass[i];//这里冲掉了原来的密度,对否??????/?
}
// Calculate SPH sum for rho:
for( k=1;k<=niac;k++)
{
i = pair_i[k];
int j = pair_j[k];
rho[i] = rho[i] + mass[j]*w[k];//这里冲掉了原来的密度,对否??????/?
rho[j] = rho[j] + mass[i]*w[k];//这里冲掉了原来的密度,对否??????/?
}
// Thirdly, calculate the normalized rho, rho=sum(rho)/sum(w)
if (nor_density) // nor_density已经定义为布尔变量；
{
for(int i=1;i<= ntotal;i++)
rho[i]=rho[i]/wi[i];//这里冲掉了原来的密度,对否??????/?
}
delete []hv;
delete []wi;
}
void SPH::int_force(int ntotal,double **(&dvx1dt),double *tdsdt,double *(&de1dt))
{
// int_force(ntotal+nvirt,indvxdt,tdsdt,du)
double *dvx1, *txx, *tyy,*tzz, *txy, *txz, *tyz, *vcc;
dvx1=new double[dim+1];txx=new double[maxn+1];tyy=new double[maxn+1];
tzz=new double[maxn+1];txy=new double[maxn+1];txz=new double[maxn+1];
tyz=new double[maxn+1];vcc=new double[maxn+1];
double hxx, hyy, hzz, hxy, hxz, hyz, h, hvcc, he, rhoij;
double temp1,temp2,temp3,temp4;
for(int i=1;i<=ntotal;i++)// Initialization of shear tensor, velocity divergence,viscous energy, internal energy, acceleration
{
txx[i] = 0.0;tyy[i] = 0.0;tzz[i] = 0.0;txy[i] = 0.0;txz[i] = 0.0;
tyz[i] = 0.0;vcc[i] = 0.0;tdsdt[i] = 0.0;de1dt[i] = 0.0;
for(int d=1;d<=dim;d++)
dvx1dt[d][i] = 0.0;
}
if (visc) // Calculate SPH sum for shear tensor Tab = va,b + vb,a - 2/3 delta_ab vc,c
{
for(int k=1;k<=niac;k++)
{
i = pair_i[k];int j = pair_j[k];
for(int d=1;d<=dim;d++)
{
dvx1[d] = vx[d][j] - vx[d][i];
}
if(dim==1)
{
hxx = 2.0*dvx1[1]*dwdx[1][k];
}
else if (dim==2)
{
hxx = 2.0*dvx1[1]*dwdx[1][k] - dvx1[2]*dwdx[2][k];
hxy = dvx1[1]*dwdx[2][k] + dvx1[2]*dwdx[1][k];
hyy = 2.0*dvx1[2]*dwdx[2][k] - dvx1[1]*dwdx[1][k];
}
else if (dim==3)
{
hxx = 2.0*dvx1[1]*dwdx[1][k] - dvx1[2]*dwdx[2][k]- dvx1[3]*dwdx[3][k];
hxy = dvx1[1]*dwdx[2][k] + dvx1[2]*dwdx[1][k];
hxz = dvx1[1]*dwdx[3][k] + dvx1[3]*dwdx[1][k];
hyy = 2.0*dvx1[2]*dwdx[2][k] - dvx1[1]*dwdx[1][k]- dvx1[3]*dwdx[3][k];
hyz = dvx1[2]*dwdx[3][k] + dvx1[3]*dwdx[2][k];
hzz = 2.0*dvx1[3]*dwdx[3][k] - dvx1[1]*dwdx[1][k]- dvx1[2]*dwdx[2][k];
}
hxx = 2.0/3.0*hxx;
hyy = 2.0/3.0*hyy;
hzz = 2.0/3.0*hzz;
if (dim==1)
{
txx[i] = txx[i] + mass[j]*hxx/rho[j];
txx[j] = txx[j] + mass[i]*hxx/rho[i];
}
else if (dim==2)
{
txx[i] = txx[i] + mass[j]*hxx/rho[j];
txx[j] = txx[j] + mass[i]*hxx/rho[i];
txy[i] = txy[i] + mass[j]*hxy/rho[j];
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -