diff --git a/src/Interaction/Models/contactForce/nonLinearCF.H b/src/Interaction/Models/contactForce/nonLinearCF.H index fbd53152..dee111ea 100644 --- a/src/Interaction/Models/contactForce/nonLinearCF.H +++ b/src/Interaction/Models/contactForce/nonLinearCF.H @@ -40,16 +40,15 @@ public: { real Yeff_ = 1000000.0; real Geff_ = 8000000.0; - real ethan_ = 0.0; - real ethat_ = 0.0; - real mu_ = 0.00001; + real ethan_= 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 ): - Yeff_(Yeff), Geff_(Geff), ethan_(etha_n),ethat_(etha_t), mu_(mu) + nonLinearProperties(real Yeff, real Geff, real etha_n, real mu ): + Yeff_(Yeff), Geff_(Geff), ethan_(etha_n), mu_(mu) {} INLINE_FUNCTION_HD @@ -78,7 +77,6 @@ protected: auto Geff = dict.getVal("Geff"); auto nu = dict.getVal("nu"); auto en = dict.getVal("en"); - auto et = dict.getVal("et"); auto mu = dict.getVal("mu"); auto nElem = Yeff.size(); @@ -97,12 +95,6 @@ protected: return false; } - if(nElem != et.size()) - { - fatalErrorInFunction<< - "sizes of Yeff("< 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(-4.0/3.0) * prop.Yeff_ * sqrt(Reff)* pow(ovrlp_n,static_cast(1.5)) - sqrt_meff_K_hertz*prop.ethan_*pow(ovrlp_n,static_cast(0.25))*vrn)*Nij; - FCt = (- static_cast(16.0/3.0) * prop.Geff_ * sqrt(Reff*ovrlp_n) ) * history.overlap_t_; + FCt = (- static_cast(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(16.0/3.0) * prop.Geff_ * sqrt(Reff*ovrlp_n); + real kt = static_cast(8.0) * prop.Geff_ * sqrt(Reff*ovrlp_n); FCt *= (ft_fric/ft); history.overlap_t_ = - (FCt/kt); } diff --git a/src/Interaction/Models/contactForce/nonLinearMod.H b/src/Interaction/Models/contactForce/nonLinearMod.H new file mode 100644 index 00000000..c471f7b0 --- /dev/null +++ b/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 +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; + + int32 numMaterial_ = 0; + + ViewType1D rho_; + + NonLinearArrayType nonlinearProperties_; + + bool readNonLinearDictionary(const dictionary& dict) + { + auto Yeff = dict.getVal("Yeff"); + auto Geff = dict.getVal("Geff"); + auto nu = dict.getVal("nu"); + auto etha_n = dict.getVal("etha_n"); + auto mu = dict.getVal("mu"); + + auto nElem = Yeff.size(); + + if(nElem != nu.size()) + { + fatalErrorInFunction<< + "sizes of Yeff("< 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& 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(3))*rho_[propId_i]; + real mj = 3*Pi/4*pow(Rj,static_cast(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(-4.0/3.0) * prop.Yeff_ * sqrt(Reff)* pow(ovrlp_n,static_cast(1.5)) - + prop.ethan_*pow(ovrlp_n,static_cast(0.5))*vrn)*Nij; + + FCt = (- static_cast(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(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__ diff --git a/src/Interaction/Models/contactForceModels.H b/src/Interaction/Models/contactForceModels.H index b359cfa2..de2117b3 100644 --- a/src/Interaction/Models/contactForceModels.H +++ b/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>; using limitedNonLinearNormalRolling = normalRolling>; using nonLimitedNonLinearNormalRolling = normalRolling>; +using limitedNonLinearModNormalRolling = normalRolling>; +using nonLimitedNonLinearModNormalRolling = normalRolling>; + } diff --git a/src/Interaction/sphereInteraction/sphereInteractions.C b/src/Interaction/sphereInteraction/sphereInteractions.C index 7fbafda7..04ef878e 100644 --- a/src/Interaction/sphereInteraction/sphereInteractions.C +++ b/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>; +