granular flow More...

#include <navierStokesAdvectionDiffusionCoupling.h>

Collaboration diagram for olb::GranularCoupling< T >:

Classes
struct	FORCE_PREFACTOR

struct	FRICTION_ANGLE

struct	RHO0

Public Types
using	parameters = meta::list<FRICTION_ANGLE,RHO0,FORCE_PREFACTOR>

Public Member Functions
template<typename CELLS , typename PARAMETERS >
void	apply (CELLS &cells, PARAMETERS &parameters) any_platform

Static Public Attributes
static constexpr OperatorScope	scope = OperatorScope::PerCellWithParameters

Detailed Description

template<typename T>
struct olb::GranularCoupling< T >

granular flow

Definition at line 465 of file navierStokesAdvectionDiffusionCoupling.h.

Member Typedef Documentation

◆ parameters

template<typename T >

using olb::GranularCoupling< T >::parameters = meta::list<FRICTION_ANGLE,RHO0,FORCE_PREFACTOR>

Definition at line 472 of file navierStokesAdvectionDiffusionCoupling.h.

Member Function Documentation

◆ apply()

template<typename T >

template<typename CELLS , typename PARAMETERS >

void olb::GranularCoupling< T >::apply	(	CELLS &	cells,
		PARAMETERS &	parameters )

inline

Definition at line 475 of file navierStokesAdvectionDiffusionCoupling.h.

   {
    using DESCRIPTOR = typename CELLS::template value_t<names::NavierStokes>::descriptor_t;
    T rho, u[DESCRIPTOR::d], pi[util::TensorVal<DESCRIPTOR>::n] { };
    cells.template get<names::NavierStokes>().computeAllMomenta(rho, u, pi);
    // Stress tensor
    T PiNeqNormSqr = pi[0]*pi[0] + 2.0*pi[1]*pi[1] + pi[2]*pi[2];
    if constexpr (util::TensorVal<DESCRIPTOR>::n == 6) {
      PiNeqNormSqr += pi[2]*pi[2] + pi[3]*pi[3] + 2*pi[4]*pi[4] +pi[5]*pi[5];
    }
    T PiNeqNorm = util::sqrt(PiNeqNormSqr);
 
    if(PiNeqNorm != T{0.}){
      auto omega = cells.template get<names::NavierStokes>().template getFieldPointer<descriptors::OMEGA>();
      PiNeqNorm /= (-rho/omega[0]/descriptors::invCs2<T,DESCRIPTOR>());
      T angle = parameters.template get<FRICTION_ANGLE>();
      T ps = ( T{0.5}*omega[0] - T{1} ) * T{1./3.} * ( pi[0] + pi[3] + pi[5] );
      T granularTau = ps/rho * util::sin(T{3.14}*angle/T{180}) / PiNeqNorm * descriptors::invCs2<T,DESCRIPTOR>() + T{0.5};
      omega[0] = T{1} / granularTau;
    }
 
    // computation of the bousinessq force
    auto force = cells.template get<names::NavierStokes>().template getFieldPointer<descriptors::FORCE>();
    auto forcePrefactor = parameters.template get<FORCE_PREFACTOR>();
    T densityDifference = (rho - parameters.template get<RHO0>()) / rho;
    for (unsigned iD = 0; iD < DESCRIPTOR::d; ++iD) {
      force[iD] = forcePrefactor[iD] * densityDifference;
    }
  }