# PhasicFlow: src/Interaction/Models/contactForce/nonLinearCF.hpp Source File

nonLinearCF.hpp
Go to the documentation of this file.
/*------------------------------- 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 __nonLinearCF_hpp__
#define __nonLinearCF_hpp__
#include "types.hpp"
namespace pFlow::cfModels
{
template<bool limited=true>
class nonLinear
{
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 en = dict.getVal<realVector>("en");
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 != en.size())
{
fatalErrorInFunction<<
"sizes of Yeff("<<nElem<<") and en("<<en.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;
}
realVector etha_n(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 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
}
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 "nonLinearLimited";
}
else
{
return "nonLinearNonLimited";
}
return "";
}
public:
TypeInfoNV(modelName());
INLINE_FUNCTION_HD
nonLinear(){}
nonLinear(
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
nonLinear(const nonLinear&) = default;
INLINE_FUNCTION_HD
nonLinear(nonLinear&&) = default;
INLINE_FUNCTION_HD
nonLinear& operator=(const nonLinear&) = default;
INLINE_FUNCTION_HD
nonLinear& operator=(nonLinear&&) = default;
INLINE_FUNCTION_HD
~nonLinear()=default;
INLINE_FUNCTION_HD
int32 numMaterial()const
{
return numMaterial_;
}
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)) -
sqrt_meff_K_hertz*prop.ethan_*pow(ovrlp_n,static_cast< real >(0.25))*vrn)*Nij;
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);
// apply friction
if(ft > ft_fric)
{
if( length(history.overlap_t_) >0.0)
{
if constexpr (limited)
{
real kt = static_cast< real >(8.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
pFlow::cfModels::nonLinear::nonLinearProperties::~nonLinearProperties
INLINE_FUNCTION_HD ~nonLinearProperties()=default
pFlow::cfModels::nonLinear::TypeInfoNV
TypeInfoNV(modelName())
pFlow::real
float real
Definition: builtinTypes.hpp:46
fatalExit
#define fatalExit
Definition: error.hpp:57
pFlow::cfModels::nonLinear::nonLinearProperties::Geff_
real Geff_
Definition: nonLinearCF.hpp:42
pFlow::cfModels::nonLinear::~nonLinear
INLINE_FUNCTION_HD ~nonLinear()=default
pFlow::cfModels::nonLinear::contactForce
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
Definition: nonLinearCF.hpp:215
types.hpp
pFlow::cfModels::nonLinear::nonLinear
INLINE_FUNCTION_HD nonLinear()
Definition: nonLinearCF.hpp:170
pFlow::uint32
unsigned int uint32
Definition: builtinTypes.hpp:59
pFlow::cfModels::nonLinear::nonLinearProperties::mu_
real mu_
Definition: nonLinearCF.hpp:44
pFlow::Vector::size
auto size() const
Definition: Vector.hpp:299
pFlow::cfModels::nonLinear::nonLinearProperties::nonLinearProperties
INLINE_FUNCTION_HD nonLinearProperties(real Yeff, real Geff, real etha_n, real mu)
Definition: nonLinearCF.hpp:50
pFlow::cfModels::nonLinear::nonLinearProperties::ethan_
real ethan_
Definition: nonLinearCF.hpp:43
pFlow::log
INLINE_FUNCTION_HD real log(real x)
Definition: math.hpp:119
pFlow::symArray< nonLinearProperties >::getN
static bool getN(uint32 nElem, uint32 &n)
Definition: symArrayHD.hpp:240
pFlow::cfModels::nonLinear::nonlinearProperties_
NonLinearArrayType nonlinearProperties_
Definition: nonLinearCF.hpp:72
dot
INLINE_FUNCTION_HD T dot(const quadruple< T > &oprnd1, const quadruple< T > &oprnd2)
pFlow::symArray::assign
bool assign(const Vector< T > src)
Definition: symArrayHD.hpp:177
fatalErrorInFunction
#define fatalErrorInFunction
Definition: error.hpp:42
length
INLINE_FUNCTION_HD T length(const triple< T > &v1)
pFlow::int32
int int32
Definition: builtinTypes.hpp:53
pFlow::pow
Vector< T, Allocator > pow(const Vector< T, Allocator > &v, T e)
Definition: VectorMath.hpp:109
pFlow::cfModels::nonLinear::rho_
ViewType1D< real > rho_
Definition: nonLinearCF.hpp:70
pFlow::cfModels::nonLinear::nonLinear
nonLinear(int32 nMaterial, const ViewType1D< real > &rho, const dictionary &dict)
Definition: nonLinearCF.hpp:172
ForAll
#define ForAll(i, container)
Definition: pFlowMacros.hpp:71
pFlow::cfModels::nonLinear::operator=
INLINE_FUNCTION_HD nonLinear & operator=(const nonLinear &)=default
pFlow::cfModels::nonLinear::contactForceStorage::overlap_t_
realx3 overlap_t_
Definition: nonLinearCF.hpp:36
pFlow::cfModels::nonLinear::modelName
static const char * modelName()
Definition: nonLinearCF.hpp:152
pFlow::cfModels::nonLinear
Definition: nonLinearCF.hpp:30
pFlow::dictionary::getVal
T getVal(const word &keyword) const
Definition: dictionary.hpp:309
pFlow::cfModels::nonLinear::nonLinearProperties
Definition: nonLinearCF.hpp:39
pFlow::ViewType1D
Kokkos::View< T *, properties... > ViewType1D
Definition: KokkosTypes.hpp:62
pFlow::cfModels::nonLinear::readNonLinearDictionary
bool readNonLinearDictionary(const dictionary &dict)
Definition: nonLinearCF.hpp:74
pFlow::Pi
const real Pi
Definition: numericConstants.hpp:32
pFlow::cfModels
Definition: linearCF.hpp:27
pFlow::cfModels::nonLinear::numMaterial
INLINE_FUNCTION_HD int32 numMaterial() const
Definition: nonLinearCF.hpp:206
pFlow::sqrt
INLINE_FUNCTION_HD real sqrt(real x)
Definition: math.hpp:148
pFlow::cfModels::nonLinear::nonLinearProperties::Yeff_
real Yeff_
Definition: nonLinearCF.hpp:41
INLINE_FUNCTION_HD
#define INLINE_FUNCTION_HD
Definition: pFlowMacros.hpp:51
pFlow::triple< real >
pFlow::Vector< real >
pFlow::cfModels::nonLinear::nonLinearProperties::operator=
INLINE_FUNCTION_HD nonLinearProperties & operator=(const nonLinearProperties &)=default
pFlow::cfModels::nonLinear::contactForceStorage
Definition: nonLinearCF.hpp:34
pFlow::cfModels::nonLinear::numMaterial_
int32 numMaterial_
Definition: nonLinearCF.hpp:68
pFlow::cfModels::nonLinear::nonLinearProperties::nonLinearProperties
INLINE_FUNCTION_HD nonLinearProperties()
Definition: nonLinearCF.hpp:47
pFlow::symArray< nonLinearProperties >
pFlow::dictionary
Definition: dictionary.hpp:38