Merge pull request #53 from PhasicFlow/conForce

nonlinear modified

src/Interaction/Models/contactForce/nonLinearCF.H

@@ -40,16 +40,15 @@ public:
 	{		
 		real 	Yeff_	= 1000000.0;
 		real 	Geff_	= 8000000.0;
-		real 	ethan_ 	= 0.0;
+		real 	ethan_= 0.0;
-		real 	ethat_ 	= 0.0;
 		real  mu_ 	= 0.00001;
 
 		INLINE_FUNCTION_HD
 		nonLinearProperties(){}
 
 		INLINE_FUNCTION_HD
-		nonLinearProperties(real Yeff, real Geff, real etha_n, real etha_t, real mu ):
+		nonLinearProperties(real Yeff, real Geff, real etha_n, real mu ):
-			Yeff_(Yeff), Geff_(Geff), ethan_(etha_n),ethat_(etha_t), mu_(mu)
+			Yeff_(Yeff), Geff_(Geff), ethan_(etha_n), mu_(mu)
 		{}		
 
 		INLINE_FUNCTION_HD
@@ -78,7 +77,6 @@ protected:
 		auto Geff = dict.getVal<realVector>("Geff");
 		auto nu = dict.getVal<realVector>("nu");
 		auto en = dict.getVal<realVector>("en");
-		auto et = dict.getVal<realVector>("et");
 		auto mu = dict.getVal<realVector>("mu");
 
 		auto nElem = Yeff.size();
@@ -97,12 +95,6 @@ protected:
 			return false;
 		}
 
-		if(nElem != et.size())
-		{
-			fatalErrorInFunction<<
-			"sizes of Yeff("<<nElem<<") and et("<<et.size()<<") do not match.\n";
-			return false;
-		}
 
 		if(nElem != mu.size())
 		{
@@ -130,26 +122,25 @@ protected:
 
 
 		realVector etha_n(nElem);
-		realVector etha_t(nElem);          
+		
 		forAll(i , en)
 		{
 			//K_hertz = 4.0/3.0*Yeff*sqrt(Reff); 
 			//-2.2664*log(en)*sqrt(meff*K_hertz)/sqrt( log(en)**2 + 10.1354);
 
-			// we take out sqrt(meff*K_hertz) here and then condier this term 
+			// we take out sqrt(meff*K_hertz) here and then consider this term 
 			// when calculating damping part. 
 			etha_n[i] = -2.2664*log(en[i])/
 					sqrt(pow(log(en[i]),2.0)+ pow(Pi,2.0));
 
 			// no damping for tangential part 
-			etha_t[i] = 0.0;
 
 		}
 
 		Vector<nonLinearProperties> prop(nElem);
 		forAll(i,Yeff)
 		{
-			prop[i] = {Yeff[i], Geff[i], etha_n[i], etha_t[i], mu[i]};
+			prop[i] = {Yeff[i], Geff[i], etha_n[i], mu[i]};
 		}
 
 		nonlinearProperties_.assign(prop);
@@ -253,16 +244,10 @@ public:
 		real K_hertz = 4.0/3.0*prop.Yeff_*sqrt(Reff);
 		real sqrt_meff_K_hertz = sqrt((mi*mj)/(mi+mj) * K_hertz);
 
-		
-		//FCn = (-prop.kn_ * ovrlp_n - sqrt_meff * prop.ethan_ * vrn)*Nij;
-		//FCt = -prop.kt_ * history.overlap_t_ - sqrt_meff * prop.ethat_*Vt;
-		
-
-
 		FCn = (static_cast<real>(-4.0/3.0) * prop.Yeff_ * sqrt(Reff)* pow(ovrlp_n,static_cast<real>(1.5)) - 
 			   sqrt_meff_K_hertz*prop.ethan_*pow(ovrlp_n,static_cast<real>(0.25))*vrn)*Nij;
 
-		FCt = (- static_cast<real>(16.0/3.0) * prop.Geff_ * sqrt(Reff*ovrlp_n) ) * history.overlap_t_;
+		FCt = (- static_cast<real>(8.0) * prop.Geff_ * sqrt(Reff*ovrlp_n) ) * history.overlap_t_;
 
 		real ft = length(FCt);
 		real ft_fric = prop.mu_ * length(FCn);
@@ -274,7 +259,7 @@ public:
 			{
 				if constexpr (limited)
 				{
-					real kt = static_cast<real>(16.0/3.0) * prop.Geff_ * sqrt(Reff*ovrlp_n);
+					real kt = static_cast<real>(8.0) * prop.Geff_ * sqrt(Reff*ovrlp_n);
 					FCt *= (ft_fric/ft);
                 	history.overlap_t_ = - (FCt/kt);
 				}

src/Interaction/Models/contactForce/nonLinearMod.H

@@ -0,0 +1,269 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+#ifndef __nonLinearModCF_H__
+#define __nonLinearModCF_H__
+
+#include "types.H"
+
+namespace pFlow::cfModels
+{
+
+template<bool limited=true>
+class nonLinearMod
+{
+public:
+
+	struct contactForceStorage
+	{
+		realx3 overlap_t_ = 0.0;
+	};
+
+	struct nonLinearProperties
+	{		
+		real 	Yeff_	= 1000000.0;
+		real 	Geff_	= 8000000.0;
+		real 	ethan_= 0.0;
+		real  mu_ 	= 0.00001;
+
+		INLINE_FUNCTION_HD
+		nonLinearProperties(){}
+
+		INLINE_FUNCTION_HD
+		nonLinearProperties(real Yeff, real Geff, real etha_n, real mu ):
+			Yeff_(Yeff), Geff_(Geff), ethan_(etha_n), mu_(mu)
+		{}		
+
+		INLINE_FUNCTION_HD
+		nonLinearProperties(const nonLinearProperties&)=default;
+
+		INLINE_FUNCTION_HD
+		nonLinearProperties& operator=(const nonLinearProperties&)=default;
+
+		INLINE_FUNCTION_HD
+		~nonLinearProperties() = default;
+	};
+
+protected:
+
+	using NonLinearArrayType = 	symArray<nonLinearProperties>;
+
+	int32 						numMaterial_ = 0;
+
+	ViewType1D<real>  		 	rho_;
+
+	NonLinearArrayType  		nonlinearProperties_;
+
+	bool readNonLinearDictionary(const dictionary& dict)
+	{
+		auto Yeff = dict.getVal<realVector>("Yeff");
+		auto Geff = dict.getVal<realVector>("Geff");
+		auto nu = dict.getVal<realVector>("nu");
+		auto etha_n = dict.getVal<realVector>("etha_n");
+		auto mu = dict.getVal<realVector>("mu");
+
+		auto nElem = Yeff.size();
+
+		if(nElem != nu.size())
+		{
+			fatalErrorInFunction<<
+			"sizes of Yeff("<<nElem<<") and nu("<<nu.size()<<") do not match.\n";
+			return false;
+		}
+
+		if(nElem != etha_n.size())
+		{
+			fatalErrorInFunction<<
+			"sizes of Yeff("<<nElem<<") and etha_n("<<etha_n.size()<<") do not match.\n";
+			return false;
+		}
+
+
+		if(nElem != mu.size())
+		{
+			fatalErrorInFunction<<
+			"sizes of Yeff("<<nElem<<") and mu("<<mu.size()<<") do not match.\n";
+			return false;
+		}
+
+		
+		// check if size of vector matchs a symetric array
+		uint32 nMat;
+		if( !NonLinearArrayType::getN(nElem, nMat) )
+		{
+			fatalErrorInFunction<<
+			"sizes of properties do not match a symetric array with size ("<<
+			numMaterial_<<"x"<<numMaterial_<<").\n";
+			return false;
+		}
+		else if( numMaterial_ != nMat)
+		{
+			fatalErrorInFunction<<
+			"size mismatch for porperties. \n";
+			return false;
+		}
+
+
+		Vector<nonLinearProperties> prop(nElem);
+		forAll(i,Yeff)
+		{
+			prop[i] = {Yeff[i], Geff[i], etha_n[i], mu[i]};
+		}
+
+		nonlinearProperties_.assign(prop);
+
+		return true;
+
+	}
+
+	static const char* modelName()
+	{
+		if constexpr (limited)
+		{
+			return "nonLinearModLimited";
+		}
+		else
+		{
+			return "nonLinearModNonLimited";
+		}
+		return "";
+	}
+
+public:
+
+	TypeNameNV(modelName());
+
+	INLINE_FUNCTION_HD
+	nonLinearMod(){}
+
+	nonLinearMod(
+		int32 nMaterial,
+		const ViewType1D<real>& rho,
+		const dictionary& dict)
+	:
+		numMaterial_(nMaterial),
+		rho_("rho",nMaterial),
+		nonlinearProperties_("nonLinearProperties",nMaterial)
+	{
+
+		Kokkos::deep_copy(rho_,rho);
+		if(!readNonLinearDictionary(dict))
+		{
+			fatalExit;
+		}
+	}
+	
+	INLINE_FUNCTION_HD
+	nonLinearMod(const nonLinearMod&) = default;
+
+	INLINE_FUNCTION_HD
+	nonLinearMod(nonLinearMod&&) = default;
+
+	INLINE_FUNCTION_HD
+	nonLinearMod& operator=(const nonLinearMod&) = default;
+
+	INLINE_FUNCTION_HD
+	nonLinearMod& operator=(nonLinearMod&&) = default;
+
+
+	INLINE_FUNCTION_HD
+	~nonLinearMod()=default;
+
+	INLINE_FUNCTION_HD
+	int32 numMaterial()const
+	{
+		return numMaterial_;
+	}
+
+	//// - Methods
+
+	INLINE_FUNCTION_HD
+	void contactForce
+	(
+		const real dt,
+		const int32 i,
+		const int32 j,
+		const int32 propId_i,
+		const int32 propId_j,
+		const real Ri,
+		const real Rj,
+		const real ovrlp_n,
+		const realx3& Vr,
+		const realx3& Nij,
+		contactForceStorage& history,
+		realx3& FCn,
+		realx3& FCt
+	)const
+	{
+
+		auto prop = nonlinearProperties_(propId_i,propId_j);
+
+
+		real vrn = dot(Vr, Nij);	
+		realx3 Vt = Vr - vrn*Nij;
+
+		history.overlap_t_ += Vt*dt;
+
+		real mi = 3*Pi/4*pow(Ri,static_cast<real>(3))*rho_[propId_i];
+		real mj = 3*Pi/4*pow(Rj,static_cast<real>(3))*rho_[propId_j];
+		real Reff = 1.0/(1/Ri + 1/Rj);
+
+		real K_hertz = 4.0/3.0*prop.Yeff_*sqrt(Reff);
+		real sqrt_meff_K_hertz = sqrt((mi*mj)/(mi+mj) * K_hertz);
+
+		FCn = (static_cast<real>(-4.0/3.0) * prop.Yeff_ * sqrt(Reff)* pow(ovrlp_n,static_cast<real>(1.5)) - 
+			   prop.ethan_*pow(ovrlp_n,static_cast<real>(0.5))*vrn)*Nij;
+
+		FCt = (- static_cast<real>(16.0/3.0) * prop.Geff_ * sqrt(Reff*ovrlp_n) ) * history.overlap_t_;
+
+		real ft = length(FCt);
+		real ft_fric = prop.mu_ * length(FCn);
+
+		// apply friction 
+		if(ft > ft_fric)
+		{
+			if( length(history.overlap_t_) >0.0)
+			{
+				if constexpr (limited)
+				{
+					real kt = static_cast<real>(16.0/3.0) * prop.Geff_ * sqrt(Reff*ovrlp_n);
+					FCt *= (ft_fric/ft);
+                	history.overlap_t_ = - (FCt/kt);
+				}
+				else
+				{
+					FCt = (FCt/ft)*ft_fric;	
+				}
+				
+			}
+			else
+			{
+				FCt = 0.0;
+			}
+		}
+
+	}
+   
+
+};
+
+} //pFlow::CFModels
+
+#endif //__nonLinearModCF_H__

src/Interaction/Models/contactForceModels.H

@@ -25,6 +25,7 @@ Licence:
 #include "linearCF.H"
 #include "nonLinearCF.H"
 #include "normalRolling.H"
+#include "nonLinearMod.H"
 
 
 namespace pFlow::cfModels 
@@ -37,6 +38,9 @@ using nonLimitedLinearNormalRolling = normalRolling<linear<false>>;
 using limitedNonLinearNormalRolling 	= normalRolling<nonLinear<true>>;
 using nonLimitedNonLinearNormalRolling 	= normalRolling<nonLinear<false>>;
 
+using limitedNonLinearModNormalRolling   = normalRolling<nonLinearMod<true>>;
+using nonLimitedNonLinearModNormalRolling  = normalRolling<nonLinearMod<false>>;
+
 }
 
 

src/Interaction/sphereInteraction/sphereInteractions.C

@@ -111,6 +111,47 @@ template class pFlow::sphereInteraction<
 	pFlow::sortedContactList>;
 
 
+// - nonLinearMod models
+	template class pFlow::sphereInteraction<
+	pFlow::cfModels::limitedNonLinearModNormalRolling,
+	pFlow::fixedGeometry,
+	pFlow::unsortedContactList>;
+
+template class pFlow::sphereInteraction<
+	pFlow::cfModels::limitedNonLinearModNormalRolling,
+	pFlow::fixedGeometry,
+	pFlow::sortedContactList>;
+
+template class pFlow::sphereInteraction<
+	pFlow::cfModels::nonLimitedNonLinearModNormalRolling,
+	pFlow::fixedGeometry,
+	pFlow::unsortedContactList>;
+
+template class pFlow::sphereInteraction<
+	pFlow::cfModels::nonLimitedNonLinearModNormalRolling,
+	pFlow::fixedGeometry,
+	pFlow::sortedContactList>;
+
+template class pFlow::sphereInteraction<
+	pFlow::cfModels::limitedNonLinearModNormalRolling,
+	pFlow::rotationAxisMotionGeometry,
+	pFlow::unsortedContactList>;
+
+template class pFlow::sphereInteraction<
+	pFlow::cfModels::limitedNonLinearModNormalRolling,
+	pFlow::rotationAxisMotionGeometry,
+	pFlow::sortedContactList>;
+
+template class pFlow::sphereInteraction<
+	pFlow::cfModels::nonLimitedNonLinearModNormalRolling,
+	pFlow::rotationAxisMotionGeometry,
+	pFlow::unsortedContactList>;
+
+template class pFlow::sphereInteraction<
+	pFlow::cfModels::nonLimitedNonLinearModNormalRolling,
+	pFlow::rotationAxisMotionGeometry,
+	pFlow::sortedContactList>;
+