(500847d4) · Commits · 郑智淋 / lammps

src/GRANULAR/fix_pour.cpp

+2 −2

Original line number	Diff line number	Diff line
		@@ -274,8 +274,8 @@ void FixPour::init()
		// check if a shear history fix exists

		fix_history = NULL;
		if (force->pair_match("gran/hooke/history") \|\|
		force->pair_match("gran/hertz/history"))
		if (force->pair_match("gran/hooke/history",1) \|\|
		force->pair_match("gran/hertz/history",1))
		for (int i = 0; i < modify->nfix; i++)
		if (strcmp(modify->fix[i]->style,"SHEAR_HISTORY") == 0)
		fix_history = (FixShearHistory *) modify->fix[i];

src/GRANULAR/fix_wall_gran.cpp

+70 −78

Original line number	Diff line number	Diff line
		@@ -31,7 +31,7 @@
		using namespace LAMMPS_NS;

		enum{XPLANE,YPLANE,ZPLANE,ZCYLINDER};
		enum{NO_HISTORY,HISTORY,HERTZIAN};
		enum{HOOKE,HOOKE_HISTORY,HERTZ_HISTORY};

		#define BIG 1.0e20

		@@ -190,25 +190,17 @@ int FixWallGran::setmask()

		void FixWallGran::init()
		{
		// set local values from Pair values

		if (force->pair == NULL)
		error->all("Fix wall/gran is incompatible with Pair style");
		double p_xkk = (double ) force->pair->extract("xkk");
		if (!p_xkk) error->all("Fix wall/gran is incompatible with Pair style");
		xkk = *p_xkk;

		// same initialization as in pair_gran_history::init_style()

		xkkt = xkk * 2.0/7.0;
		double gammas = 0.5*gamman;
		dt = update->dt;

		// set pairstyle from granular pair style

		if (force->pair_match("gran/no_history")) pairstyle = NO_HISTORY;
		else if (force->pair_match("gran/history")) pairstyle = HISTORY;
		else if (force->pair_match("gran/hertzian")) pairstyle = HERTZIAN;
		if (force->pair_match("gran/hooke",1))
		pairstyle = HOOKE;
		else if (force->pair_match("gran/hooke/history",1))
		pairstyle = HOOKE_HISTORY;
		else if (force->pair_match("gran/hertz/history",1))
		pairstyle = HERTZ_HISTORY;
		else error->all("Fix wall/gran is incompatible with Pair style");
		}

		/* ---------------------------------------------------------------------- */
		@@ -295,20 +287,20 @@ void FixWallGran::post_force(int vflag)
		rsq = dxdx + dydy + dz*dz;

		if (rsq > radius[i]*radius[i]) {
		if (pairstyle != NO_HISTORY) {
		if (pairstyle != HOOKE) {
		shear[i][0] = 0.0;
		shear[i][1] = 0.0;
		shear[i][2] = 0.0;
		}
		} else {
		if (pairstyle == NO_HISTORY)
		no_history(rsq,dx,dy,dz,vwall,v[i],f[i],omega[i],torque[i],
		if (pairstyle == HOOKE)
		hooke(rsq,dx,dy,dz,vwall,v[i],f[i],omega[i],torque[i],
		radius[i],rmass[i]);
		else if (pairstyle == HISTORY)
		history(rsq,dx,dy,dz,vwall,v[i],f[i],omega[i],torque[i],
		else if (pairstyle == HOOKE_HISTORY)
		hooke_history(rsq,dx,dy,dz,vwall,v[i],f[i],omega[i],torque[i],
		radius[i],rmass[i],shear[i]);
		else if (pairstyle == HERTZIAN)
		hertzian(rsq,dx,dy,dz,vwall,v[i],f[i],omega[i],torque[i],
		else if (pairstyle == HERTZ_HISTORY)
		hertz_history(rsq,dx,dy,dz,vwall,v[i],f[i],omega[i],torque[i],
		radius[i],rmass[i],shear[i]);
		}
		}
		@@ -317,13 +309,13 @@ void FixWallGran::post_force(int vflag)

		/* ---------------------------------------------------------------------- */

		void FixWallGran::no_history(double rsq, double dx, double dy, double dz,
		void FixWallGran::hooke(double rsq, double dx, double dy, double dz,
		double vwall, double v,
		double f, double omega, double *torque,
		double radius, double mass)
		{
		double r,vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
		double wr1,wr2,wr3,xmeff,damp,ccel,vtr1,vtr2,vtr3,vrel;
		double wr1,wr2,wr3,meff,damp,ccel,vtr1,vtr2,vtr3,vrel;
		double fn,fs,ft,fs1,fs2,fs3,fx,fy,fz,tor1,tor2,tor3,rinv,rsqinv;

		r = sqrt(rsq);
		@@ -357,9 +349,9 @@ void FixWallGran::no_history(double rsq, double dx, double dy, double dz,

		// normal forces = Hookian contact + normal velocity damping

		xmeff = mass;
		damp = xmeffgammanvnnr*rsqinv;
		ccel = xkk(radius-r)rinv - damp;
		meff = mass;
		damp = meffgammanvnnr*rsqinv;
		ccel = kn(radius-r)rinv - damp;

		// relative velocities

		@@ -372,7 +364,7 @@ void FixWallGran::no_history(double rsq, double dx, double dy, double dz,
		// force normalization

		fn = xmu * fabs(ccel*r);
		fs = xmeffgammasvrel;
		fs = meffgammatvrel;
		if (vrel != 0.0) ft = MIN(fn,fs) / vrel;
		else ft = 0.0;

		@@ -402,13 +394,13 @@ void FixWallGran::no_history(double rsq, double dx, double dy, double dz,

		/* ---------------------------------------------------------------------- */

		void FixWallGran::history(double rsq, double dx, double dy, double dz,
		void FixWallGran::hooke_history(double rsq, double dx, double dy, double dz,
		double vwall, double v,
		double f, double omega, double *torque,
		double radius, double mass, double *shear)
		{
		double r,vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
		double wr1,wr2,wr3,xmeff,damp,ccel,vtr1,vtr2,vtr3,vrel;
		double wr1,wr2,wr3,meff,damp,ccel,vtr1,vtr2,vtr3,vrel;
		double fn,fs,fs1,fs2,fs3,fx,fy,fz,tor1,tor2,tor3;
		double shrmag,rsht,rinv,rsqinv;

		@@ -443,9 +435,9 @@ void FixWallGran::history(double rsq, double dx, double dy, double dz,

		// normal forces = Hookian contact + normal velocity damping

		xmeff = mass;
		damp = xmeffgammanvnnr*rsqinv;
		ccel = xkk(radius-r)rinv - damp;
		meff = mass;
		damp = meffgammanvnnr*rsqinv;
		ccel = kn(radius-r)rinv - damp;

		// relative velocities

		@@ -472,9 +464,9 @@ void FixWallGran::history(double rsq, double dx, double dy, double dz,

		// tangential forces = shear + tangential velocity damping

		fs1 = - (xkktshear[0] + xmeffgammas*vtr1);
		fs2 = - (xkktshear[1] + xmeffgammas*vtr2);
		fs3 = - (xkktshear[2] + xmeffgammas*vtr3);
		fs1 = - (ktshear[0] + meffgammat*vtr1);
		fs2 = - (ktshear[1] + meffgammat*vtr2);
		fs3 = - (ktshear[2] + meffgammat*vtr3);

		// rescale frictional displacements and forces if needed

		@@ -483,15 +475,15 @@ void FixWallGran::history(double rsq, double dx, double dy, double dz,

		if (fs > fn) {
		if (shrmag != 0.0) {
		shear[0] = (fn/fs) * (shear[0] + xmeffgammasvtr1/xkkt) -
		xmeffgammasvtr1/xkkt;
		shear[1] = (fn/fs) * (shear[1] + xmeffgammasvtr2/xkkt) -
		xmeffgammasvtr2/xkkt;
		shear[2] = (fn/fs) * (shear[2] + xmeffgammasvtr3/xkkt) -
		xmeffgammasvtr3/xkkt;
		fs1 = fs1 * fn / fs ;
		fs2 = fs2 * fn / fs;
		fs3 = fs3 * fn / fs;
		shear[0] = (fn/fs) * (shear[0] + meffgammatvtr1/kt) -
		meffgammatvtr1/kt;
		shear[1] = (fn/fs) * (shear[1] + meffgammatvtr2/kt) -
		meffgammatvtr2/kt;
		shear[2] = (fn/fs) * (shear[2] + meffgammatvtr3/kt) -
		meffgammatvtr3/kt;
		fs1 *= fn/fs ;
		fs2 *= fn/fs;
		fs3 *= fn/fs;
		} else fs1 = fs2 = fs3 = 0.0;
		}

		@@ -515,15 +507,15 @@ void FixWallGran::history(double rsq, double dx, double dy, double dz,

		/* ---------------------------------------------------------------------- */

		void FixWallGran::hertzian(double rsq, double dx, double dy, double dz,
		void FixWallGran::hertz_history(double rsq, double dx, double dy, double dz,
		double vwall, double v,
		double f, double omega, double *torque,
		double radius, double mass, double *shear)
		{
		double r,vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
		double wr1,wr2,wr3,xmeff,damp,ccel,vtr1,vtr2,vtr3,vrel;
		double wr1,wr2,wr3,meff,damp,ccel,vtr1,vtr2,vtr3,vrel;
		double fn,fs,fs1,fs2,fs3,fx,fy,fz,tor1,tor2,tor3;
		double shrmag,rsht,rhertz,rinv,rsqinv;
		double shrmag,rsht,polyhertz,rinv,rsqinv;

		r = sqrt(rsq);
		rinv = 1.0/r;
		@@ -556,11 +548,11 @@ void FixWallGran::hertzian(double rsq, double dx, double dy, double dz,

		// normal forces = Hertzian contact + normal velocity damping

		xmeff = mass;
		damp = xmeffgammanvnnr*rsqinv;
		ccel = xkk(radius-r)rinv - damp;
		rhertz = sqrt(radius - r);
		ccel = rhertz * ccel;
		meff = mass;
		damp = meffgammanvnnr*rsqinv;
		ccel = kn(radius-r)rinv - damp;
		polyhertz = sqrt((radius-r)*radius);
		ccel *= polyhertz;

		// relative velocities

		@@ -587,9 +579,9 @@ void FixWallGran::hertzian(double rsq, double dx, double dy, double dz,

		// tangential forces = shear + tangential velocity damping

		fs1 = -rhertz * (xkktshear[0] + xmeffgammas*vtr1);
		fs2 = -rhertz * (xkktshear[1] + xmeffgammas*vtr2);
		fs3 = -rhertz * (xkktshear[2] + xmeffgammas*vtr3);
		fs1 = -polyhertz * (ktshear[0] + meffgammat*vtr1);
		fs2 = -polyhertz * (ktshear[1] + meffgammat*vtr2);
		fs3 = -polyhertz * (ktshear[2] + meffgammat*vtr3);

		// rescale frictional displacements and forces if needed

		@@ -598,15 +590,15 @@ void FixWallGran::hertzian(double rsq, double dx, double dy, double dz,

		if (fs > fn) {
		if (shrmag != 0.0) {
		shear[0] = (fn/fs) * (shear[0] + xmeffgammasvtr1/xkkt) -
		xmeffgammasvtr1/xkkt;
		shear[1] = (fn/fs) * (shear[1] + xmeffgammasvtr2/xkkt) -
		xmeffgammasvtr2/xkkt;
		shear[2] = (fn/fs) * (shear[2] + xmeffgammasvtr3/xkkt) -
		xmeffgammasvtr3/xkkt;
		fs1 = fs1 * fn / fs ;
		fs2 = fs2 * fn / fs;
		fs3 = fs3 * fn / fs;
		shear[0] = (fn/fs) * (shear[0] + meffgammatvtr1/kt) -
		meffgammatvtr1/kt;
		shear[1] = (fn/fs) * (shear[1] + meffgammatvtr2/kt) -
		meffgammatvtr2/kt;
		shear[2] = (fn/fs) * (shear[2] + meffgammatvtr3/kt) -
		meffgammatvtr3/kt;
		fs1 *= fn/fs ;
		fs2 *= fn/fs;
		fs3 *= fn/fs;
		} else fs1 = fs2 = fs3 = 0.0;
		}

src/GRANULAR/fix_wall_gran.h

+9 −9

Original line number	Diff line number	Diff line
		@@ -40,7 +40,7 @@ class FixWallGran : public Fix {

		private:
		int wallstyle,pairstyle,wiggle,wshear,axis;
		double xkk,xkkt,gamman,gammas,xmu;
		double kn,kt,gamman,gammat,xmu;
		double lo,hi,cylradius;
		double dt;
		double amplitude,period,omega,time_origin,vshear;
		@@ -48,12 +48,12 @@ class FixWallGran : public Fix {
		int *touch;
		double **shear;

		void no_history(double, double, double, double, double *,
		void hooke(double, double, double, double, double *,
		double , double , double , double , double, double);
		void history(double, double, double, double, double *,
		void hooke_history(double, double, double, double, double *,
		double , double , double , double , double, double,
		double *);
		void hertzian(double, double, double, double, double *,
		void hertz_history(double, double, double, double, double *,
		double , double , double , double , double, double,
		double *);
		};

src/GRANULAR/pair_gran_hertz_history.cpp

+51 −23

Original line number	Diff line number	Diff line
		@@ -16,12 +16,14 @@
		------------------------------------------------------------------------- */

		#include "math.h"
		#include "stdlib.h"
		#include "stdio.h"
		#include "string.h"
		#include "pair_gran_hertz_history.h"
		#include "atom.h"
		#include "force.h"
		#include "neigh_list.h"
		#include "error.h"

		using namespace LAMMPS_NS;

		@@ -47,9 +49,9 @@ void PairGranHertzHistory::compute(int eflag, int vflag)
		double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
		double wr1,wr2,wr3;
		double vtr1,vtr2,vtr3,vrel;
		double xmeff,damp,ccel,tor1,tor2,tor3;
		double meff,damp,ccel,tor1,tor2,tor3;
		double fn,fs,fs1,fs2,fs3;
		double shrmag,rsht,rhertz;
		double shrmag,rsht,polyhertz;
		int ilist,jlist,numneigh,*firstneigh;
		int touch,*firsttouch;
		double shear,allshear,**firstshear;
		@@ -139,13 +141,13 @@ void PairGranHertzHistory::compute(int eflag, int vflag)

		// normal force = Hertzian contact + normal velocity damping

		xmeff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]);
		if (mask[i] & freeze_group_bit) xmeff = rmass[j];
		if (mask[j] & freeze_group_bit) xmeff = rmass[i];
		damp = xmeffgammanvnnr*rsqinv;
		ccel = xkk(radsum-r)rinv - damp;
		rhertz = sqrt(radsum - r);
		ccel *= rhertz;
		meff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]);
		if (mask[i] & freeze_group_bit) meff = rmass[j];
		if (mask[j] & freeze_group_bit) meff = rmass[i];
		damp = meffgammanvnnr*rsqinv;
		ccel = kn(radsum-r)rinv - damp;
		polyhertz = sqrt((radsum-r)radiradj / radsum);
		ccel *= polyhertz;

		// relative velocities

		@@ -175,9 +177,9 @@ void PairGranHertzHistory::compute(int eflag, int vflag)

		// tangential forces = shear + tangential velocity damping

		fs1 = -rhertz * (xkktshear[0] + xmeffgammas*vtr1);
		fs2 = -rhertz * (xkktshear[1] + xmeffgammas*vtr2);
		fs3 = -rhertz * (xkktshear[2] + xmeffgammas*vtr3);
		fs1 = -polyhertz * (ktshear[0] + meffgammat*vtr1);
		fs2 = -polyhertz * (ktshear[1] + meffgammat*vtr2);
		fs3 = -polyhertz * (ktshear[2] + meffgammat*vtr3);

		// rescale frictional displacements and forces if needed

		@@ -186,20 +188,16 @@ void PairGranHertzHistory::compute(int eflag, int vflag)

		if (fs > fn) {
		if (shrmag != 0.0) {
		shear[0] = (fn/fs) * (shear[0] + xmeffgammasvtr1/xkkt) -
		xmeffgammasvtr1/xkkt;
		shear[1] = (fn/fs) * (shear[1] + xmeffgammasvtr2/xkkt) -
		xmeffgammasvtr2/xkkt;
		shear[2] = (fn/fs) * (shear[2] + xmeffgammasvtr3/xkkt) -
		xmeffgammasvtr3/xkkt;
		shear[0] = (fn/fs) * (shear[0] + meffgammatvtr1/kt) -
		meffgammatvtr1/kt;
		shear[1] = (fn/fs) * (shear[1] + meffgammatvtr2/kt) -
		meffgammatvtr2/kt;
		shear[2] = (fn/fs) * (shear[2] + meffgammatvtr3/kt) -
		meffgammatvtr3/kt;
		fs1 *= fn/fs;
		fs2 *= fn/fs;
		fs3 *= fn/fs;
		} else {
		fs1 = 0.0;
		fs2 = 0.0;
		fs3 = 0.0;
		}
		} else fs1 = fs2 = fs3 = 0.0;
		}

		// forces & torques
		@@ -233,3 +231,33 @@ void PairGranHertzHistory::compute(int eflag, int vflag)
		}
		}
		}

		/* ----------------------------------------------------------------------
		global settings
		------------------------------------------------------------------------- */

		void PairGranHertzHistory::settings(int narg, char **arg)
		{
		if (narg != 6) error->all("Illegal pair_style command");

		kn = atof(arg[0]);
		if (strcmp(arg[1],"NULL") == 0) kt = kn * 2.0/7.0;
		else kt = atof(arg[1]);

		gamman = atof(arg[2]);
		if (strcmp(arg[3],"NULL") == 0) gammat = 0.5 * gamman;
		else gammat = atof(arg[3]);

		xmu = atof(arg[4]);
		dampflag = atoi(arg[5]);
		if (dampflag == 0) gammat = 0.0;

		if (kn < 0.0 \|\| kt < 0.0 \|\| gamman < 0.0 \|\| gammat < 0.0 \|\|
		xmu < 0.0 \|\| xmu > 1.0 \|\| dampflag < 0 \|\| dampflag > 1)
		error->all("Illegal pair_style command");

		// convert Kn and Kt from pressure units to force/distance^2

		kn /= force->nktv2p;
		kt /= force->nktv2p;
		}

src/GRANULAR/pair_gran_hertz_history.h

+1 −0

Original line number	Diff line number	Diff line
		@@ -22,6 +22,7 @@ class PairGranHertzHistory : public PairGranHookeHistory {
		public:
		PairGranHertzHistory(class LAMMPS *);
		void compute(int, int);
		void settings(int, char **);
		};

		}

Admin message