From 967bb753aa72270a3c0e028dfdf4b8d195aeb6dc Mon Sep 17 00:00:00 2001 From: Hamidreza Norouzi Date: Mon, 20 Jan 2025 15:37:48 +0330 Subject: [PATCH] adding non-linear contact force model for grains --- .../Models/contactForce/cGNonLinearCF.hpp | 307 ++++++++++++++++++ .../contactForce/cGRelativeLinearCF.hpp | 60 +--- .../Models/grainContactForceModels.hpp | 4 + .../grainInteractionsNonLinearModels.cpp | 6 + 4 files changed, 334 insertions(+), 43 deletions(-) create mode 100644 src/Interaction/Models/contactForce/cGNonLinearCF.hpp diff --git a/src/Interaction/Models/contactForce/cGNonLinearCF.hpp b/src/Interaction/Models/contactForce/cGNonLinearCF.hpp new file mode 100644 index 00000000..a167f8fa --- /dev/null +++ b/src/Interaction/Models/contactForce/cGNonLinearCF.hpp @@ -0,0 +1,307 @@ +/*------------------------------- 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 __cGNonLinearCF_hpp__ +#define __cGNonLinearCF_hpp__ + +#include "types.hpp" + +namespace pFlow::cfModels +{ + +template +class cGNonLinear +{ +public: + + struct contactForceStorage + { + realx3 overlap_t_ = 0.0; + }; + + struct cGNonLinearProperties + { + real Yeff_ = 1000000.0; + real Geff_ = 8000000.0; + real en_ = 1.0; + real mu_ = 0.00001; + + INLINE_FUNCTION_HD + cGNonLinearProperties(){} + + INLINE_FUNCTION_HD + cGNonLinearProperties(real Yeff, real Geff, real en, real mu ): + Yeff_(Yeff), Geff_(Geff), en_(en), mu_(mu) + {} + + INLINE_FUNCTION_HD + cGNonLinearProperties(const cGNonLinearProperties&)=default; + + INLINE_FUNCTION_HD + cGNonLinearProperties& operator=(const cGNonLinearProperties&)=default; + + INLINE_FUNCTION_HD + ~cGNonLinearProperties() = default; + }; + +protected: + + using cGNonLinearArrayType = symArray; + + int32 numMaterial_ = 0; + + ViewType1D rho_; + + int32 addDissipationModel_; + + cGNonLinearArrayType nonlinearProperties_; + + bool readNonLinearDictionary(const dictionary& dict) + { + auto Yeff = dict.getVal("Yeff"); + auto Geff = dict.getVal("Geff"); + auto nu = dict.getVal("nu"); + auto en = dict.getVal("en"); + auto mu = dict.getVal("mu"); + + auto nElem = Yeff.size(); + + if(nElem != nu.size()) + { + fatalErrorInFunction<< + "sizes of Yeff("< prop("prop",nElem); + ForAll(i,Yeff) + { + prop[i] = {Yeff[i], Geff[i], en[i], mu[i]}; + } + + nonlinearProperties_.assign(prop); + + auto adm = dict.getVal("additionalDissipationModel"); + if(adm == "none") + { + addDissipationModel_ = 1; + } + else if(adm == "LU") + { + addDissipationModel_ = 2; + } + else if (adm == "GB") + { + addDissipationModel_ = 3; + } + else + { + addDissipationModel_ = 1; + } + + return true; + + } + + static const char* modelName() + { + if constexpr (limited) + { + return "cGNonLinearLimited"; + } + else + { + return "cGNonLinearNonLimited"; + } + return ""; + } + +public: + + TypeInfoNV(modelName()); + + INLINE_FUNCTION_HD + cGNonLinear(){} + + cGNonLinear( + int32 nMaterial, + const ViewType1D& rho, + const dictionary& dict) + : + numMaterial_(nMaterial), + rho_("rho",nMaterial), + nonlinearProperties_("cGNonLinearProperties",nMaterial) + { + + Kokkos::deep_copy(rho_,rho); + if(!readNonLinearDictionary(dict)) + { + fatalExit; + } + } + + INLINE_FUNCTION_HD + cGNonLinear(const cGNonLinear&) = default; + + INLINE_FUNCTION_HD + cGNonLinear(cGNonLinear&&) = default; + + INLINE_FUNCTION_HD + cGNonLinear& operator=(const cGNonLinear&) = default; + + INLINE_FUNCTION_HD + cGNonLinear& operator=(cGNonLinear&&) = default; + + + INLINE_FUNCTION_HD + ~cGNonLinear()=default; + + INLINE_FUNCTION_HD + int32 numMaterial()const + { + return numMaterial_; + } + + //// - Methods + + INLINE_FUNCTION_HD + void contactForce + ( + const real dt, + const uint32 i, + const uint32 j, + const uint32 propId_i, + const uint32 propId_j, + const real Ri, + const real Rj, + const real cGFi, + const real cGFj, + const real ovrlp_n, + const realx3& Vr, + const realx3& Nij, + contactForceStorage& history, + realx3& FCn, + realx3& FCt + )const + { + + const auto prop = nonlinearProperties_(propId_i,propId_j); + + const real f = 2.0/( 1.0/cGFi + 1.0/cGFj ); + + 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); + + real en = prop.en_; + if (addDissipationModel_==2) + { + en = sqrt(1+((pow(prop.en_,2)-1)*f)); + } + else if (addDissipationModel_==3) + { + + en = exp((pow(f,1.5)*log(prop.en_)*sqrt( (1-((pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(Pi,2))))/(1-(pow(f,3)*(pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(Pi,2)))) ) )); + } + + real Kn = static_cast(4.0/3.0) * prop.Yeff_ * sqrt(Reff*ovrlp_n); + + real ethan_ = -2.0*log(en)*sqrt(Kn)/ + sqrt(pow(log(en),2.0)+ pow(Pi,2.0)); + + FCn = ( - Kn*ovrlp_n - + sqrt_meff_K_hertz*ethan_*pow(ovrlp_n,static_cast(0.25))*vrn)*Nij; + + FCt = (f*f)*(- static_cast(8.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(8.0) * prop.Geff_ * sqrt(Reff*ovrlp_n); + FCt *= (ft_fric/ft); + history.overlap_t_ = - FCt/(f*f*kt); + } + else + { + FCt = (FCt/ft)*ft_fric; + } + + } + else + { + FCt = 0.0; + } + } + + } + + +}; + +} //pFlow::CFModels + +#endif diff --git a/src/Interaction/Models/contactForce/cGRelativeLinearCF.hpp b/src/Interaction/Models/contactForce/cGRelativeLinearCF.hpp index 48a77884..d542e70c 100755 --- a/src/Interaction/Models/contactForce/cGRelativeLinearCF.hpp +++ b/src/Interaction/Models/contactForce/cGRelativeLinearCF.hpp @@ -26,11 +26,6 @@ Licence: - - - - - namespace pFlow::cfModels { @@ -49,16 +44,15 @@ public: { real kn_ = 1000.0; real kt_ = 800.0; - real en_ = 0.0; - real ethat_ = 0.0; - real mu_ = 0.00001; + real en_ = 1.0; + real mu_ = 0.00001; INLINE_FUNCTION_HD linearProperties(){} INLINE_FUNCTION_HD - linearProperties(real kn, real kt, real en, real etha_t, real mu ): - kn_(kn), kt_(kt), en_(en),ethat_(etha_t), mu_(mu) + linearProperties(real kn, real kt, real en, real mu ): + kn_(kn), kt_(kt), en_(en), mu_(mu) {} INLINE_FUNCTION_HD @@ -93,7 +87,6 @@ protected: auto kn = dict.getVal("kn"); auto kt = dict.getVal("kt"); auto en = dict.getVal("en"); - auto et = dict.getVal("et"); auto mu = dict.getVal("mu"); @@ -114,12 +107,6 @@ protected: return false; } - if(nElem != et.size()) - { - fatalErrorInFunction<< - "sizes of kn("< prop("prop", nElem); ForAll(i,kn) { - prop[i] = {kn[i], kt[i], en[i], etha_t[i], mu[i] }; + prop[i] = {kn[i], kt[i], en[i], mu[i] }; } linearProperties_.assign(prop); - auto adm = dict.getVal("additionalDissipationModel"); - + auto adm = dict.getVal("additionalDissipationModel"); if(adm == "none") { addDissipationModel_ = 1; @@ -287,22 +260,23 @@ public: real sqrt_meff = sqrt((mi*mj)/(mi+mj)); - - if (addDissipationModel_==2) + real en = prop.en_; + if (addDissipationModel_==2) { - prop.en_ = sqrt(1+((pow(prop.en_,2)-1)*f_)); + en = sqrt(1+((pow(prop.en_,2)-1)*f_)); } - else if (addDissipationModel_==3) + else if (addDissipationModel_==3) { - auto pie =3.14; - prop.en_ = exp((pow(f_,1.5)*log(prop.en_)*sqrt( (1-((pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(pie,2))))/(1-(pow(f_,3)*(pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(pie,2)))) ) )); + + en = exp((pow(f_,1.5)*log(prop.en_)*sqrt( (1-((pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(Pi,2))))/(1-(pow(f_,3)*(pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(Pi,2)))) ) )); } - real ethan_ = -2.0*log(prop.en_)*sqrt(prop.kn_)/ - sqrt(pow(log(prop.en_),2.0)+ pow(Pi,2.0)); + + real ethan_ = -2.0*log(en)*sqrt(prop.kn_)/ + sqrt(pow(log(en),2.0)+ pow(Pi,2.0)); - FCn = ( -pow(f_,1.0)*prop.kn_ * ovrlp_n - sqrt_meff * pow(f_,0.5) * ethan_ * vrn)*Nij; - FCt = ( -pow(f_,1.0)*prop.kt_ * history.overlap_t_ - sqrt_meff * pow(f_,0.5) * prop.ethat_*Vt); + FCn = ( -f_*prop.kn_ * ovrlp_n - sqrt_meff * pow(f_,0.5) * ethan_ * vrn)*Nij; + FCt = ( -f_*prop.kt_ * history.overlap_t_); diff --git a/src/Interaction/Models/grainContactForceModels.hpp b/src/Interaction/Models/grainContactForceModels.hpp index 908117f5..1b4d4cbe 100755 --- a/src/Interaction/Models/grainContactForceModels.hpp +++ b/src/Interaction/Models/grainContactForceModels.hpp @@ -23,6 +23,7 @@ Licence: #include "cGAbsoluteLinearCF.hpp" #include "cGRelativeLinearCF.hpp" +#include "cGNonLinearCF.hpp" #include "grainRolling.hpp" @@ -38,6 +39,9 @@ using nonLimitedCGAbsoluteLinearGrainRolling = grainRolling>; using nonLimitedCGRelativeLinearGrainRolling = grainRolling>; +using limitedCGNonLinearGrainRolling = grainRolling>; +using nonLimitedCGNonLinearGrainRolling = grainRolling>; + } diff --git a/src/Interaction/grainInteraction/grainInteractionsNonLinearModels.cpp b/src/Interaction/grainInteraction/grainInteractionsNonLinearModels.cpp index 462b5eb3..42f0b469 100755 --- a/src/Interaction/grainInteraction/grainInteractionsNonLinearModels.cpp +++ b/src/Interaction/grainInteraction/grainInteractionsNonLinearModels.cpp @@ -41,3 +41,9 @@ Licence: createBoundaryGrainInteraction(ForceModel, GeomModel) +createInteraction(pFlow::cfModels::limitedCGNonLinearGrainRolling, pFlow::rotationAxisMotionGeometry); + +createInteraction(pFlow::cfModels::nonLimitedCGNonLinearGrainRolling, pFlow::rotationAxisMotionGeometry); + +createInteraction(pFlow::cfModels::limitedCGNonLinearGrainRolling, pFlow::stationaryGeometry); +createInteraction(pFlow::cfModels::nonLimitedCGNonLinearGrainRolling, pFlow::stationaryGeometry); \ No newline at end of file