颗粒运动源程序.DOCVIP

PROGRAM PARTICLES_MOVEMENT !二维颗粒运动示例程序 IMPLICIT NONE INTEGER,PARAMETER::NX=20,KPMAX=10000 !最大颗粒个数，最大计算步数 REAL,PARAMETER::PI=3.1415926535,G=0 !9.80665 REAL::DP,DT,RHOP,Ct,Cn,XMIN,XMAX,YMIN,YMAX,AMU,TAO,R,UT,UG,VG INTEGER::IPMAX,IP,KP,KPEND REAL::UP(NX,0:KPMAX),VP(NX,0:KPMAX),XP(NX,0:KPMAX),YP(NX,0:KPMAX) REAL::AK1,AK2,AK3,AK4 DP=100.E-6 !颗粒粒径 DT=1.E-2 !计算时间步长 RHOP=1.E3 !颗粒密度 Ct=1. !切向碰撞系数 Cn=1. !法向碰撞系数 IPMAX=1 !最大颗粒数为1，只计算一个颗粒 XMIN=-1. !计算区域最小X坐标 XMAX=1. !计算区域最大X坐标 YMIN=-1. !计算区域最小Y坐标 YMAX=1. !计算区域最大Y坐标 AMU=1.5E-5 !气体的动力粘性 DO IP=1,IPMAX XP(IP,0)=0.1 !初始的颗粒X坐标 YP(IP,0)=0. !初始的颗粒Y坐标 ENDDO TAO=RHOP*DP*DP/(18*AMU) !计算松弛时间 KPEND=0 DO IP=1,IPMAX !IP是颗粒的编号 UP(IP,0)=0. !设定颗粒的初始U速度 VP(IP,0)=0. !设定颗粒的初始V速度 DO KP=0,KPMAX-1 !KP是每个颗粒计算步数的编号 R=SQRT(XP(IP,KP)**2+YP(IP,KP)**2) UT=1.*R UG=UT*(-YP(IP,KP)/R) !设定一个刚体旋转气流流场的U速度 VG=UT*(XP(IP,KP)/R) !设定一个刚体旋转气流流场的V速度 !使用上一时层颗粒的速度和位置，用4阶R-K法推算本时层颗粒的速度和位置 AK1=DT*(UG-UP(IP,KP))/TAO AK2=DT*(UG-UP(IP,KP)-.5*AK1)/TAO AK3=DT*(UG-UP(IP,KP)-.5*AK2)/TAO AK4=DT*(UG-UP(IP,KP)-AK3)/TAO UP(IP,KP+1)=UP(IP,KP)+(AK1+AK2*2+AK3*2+AK4)/6 !4阶R-K法算颗粒的U速度 AK1=DT*((VG-VP(IP,KP))/TAO-G) AK2=DT*((VG-VP(IP,KP)-.5*AK1)/TAO-G) AK3=DT*((VG-VP(IP,KP)-.5*AK2)/TAO-G) AK4=DT*((VG-VP(IP,KP)-AK3)/TAO-G) VP(IP,KP+1)=VP(IP,KP)+(AK1+AK2*2+AK3*2+AK4)/6 !4阶R-K法算颗粒的V速度 XP(IP,KP+1)=XP(IP,KP)+(UP(IP,KP+1)+UP(IP,KP))*DT*.5 !用U速度积分X坐标值 YP(IP,KP+1)=YP(IP,KP)+(VP(IP,KP+1)+VP(IP,KP))*DT*.5 !用V速度积分Y坐标值 IF(XP(IP,KP).GT.XMAX.OR.XP(IP,KP).LT.XMIN) THEN UP(IP,KP+1)=-Cn*UP(IP,KP-1) VP(IP,KP+1)=Ct*VP(IP,KP-1) XP(IP,KP+1)=XP(IP,KP-1) YP(IP,KP+1)=YP(IP,KP-1) !颗粒碰壁回弹的模型 ENDIF IF(YP(IP,KP).GT.YMAX.OR.YP(IP,KP).LT.YMIN) THEN UP(IP,KP+1)=Ct*UP(IP,KP-1) VP(IP,KP+1)=-Cn*VP(IP,KP-1) XP(IP,KP+1)=XP(IP,KP-1) YP(IP,KP+1)=YP(IP