VANS-ADE coupling. More...

#include <navierStokesAdvectionDiffusionCoupling.h>

Collaboration diagram for olb::VANSADECoupling< T >:

Classes
struct	CONV_DENS

struct	CONV_FORCE

struct	CONV_MASS

struct	CONV_VEL

struct	DT

struct	EARTH_ACC

struct	PART_DENS

struct	PARTICLE_DIAMETER

struct	VISCOSITY

Public Types
using	parameters = meta::list<PARTICLE_DIAMETER,VISCOSITY,CONV_VEL,CONV_DENS,PART_DENS,DT,CONV_FORCE,CONV_MASS,EARTH_ACC>

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::VANSADECoupling< T >

VANS-ADE coupling.

Definition at line 353 of file navierStokesAdvectionDiffusionCoupling.h.

Member Typedef Documentation

◆ parameters

template<typename T >

using olb::VANSADECoupling< T >::parameters = meta::list<PARTICLE_DIAMETER,VISCOSITY,CONV_VEL,CONV_DENS,PART_DENS,DT,CONV_FORCE,CONV_MASS,EARTH_ACC>

Definition at line 366 of file navierStokesAdvectionDiffusionCoupling.h.

Member Function Documentation

◆ apply()

template<typename T >

template<typename CELLS , typename PARAMETERS >

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

inline

Definition at line 369 of file navierStokesAdvectionDiffusionCoupling.h.

   {
    T particleDiam = parameters.template get<PARTICLE_DIAMETER>();
    T visc = parameters.template get<VISCOSITY>();
    T convVel = parameters.template get<CONV_VEL>();
    T convDens = parameters.template get<CONV_DENS>();
    T convForce = parameters.template get<CONV_FORCE>();
    T convMass = parameters.template get<CONV_MASS>();
    T partDens = parameters.template get<PART_DENS>();
    T dt = parameters.template get<DT>();
    auto earthAcc = parameters.template get<EARTH_ACC>();
    using DESCRIPTOR = typename CELLS::template value_t<names::NavierStokes>::descriptor_t;
    auto& cellNSE = cells.template get<names::NavierStokes>();
    auto& cellADE = cells.template get<names::Concentration0>();
 
    auto u_p = cellADE.template getField<descriptors::VELOCITY>();
    auto u_pXp = cellADE.neighbor({1,0,0}).template getField<descriptors::VELOCITY2>();
    auto u_pXm = cellADE.neighbor({-1,0,0}).template getField<descriptors::VELOCITY2>();
    auto u_pYp = cellADE.neighbor({0,1,0}).template getField<descriptors::VELOCITY2>();
    auto u_pYm = cellADE.neighbor({0,-1,0}).template getField<descriptors::VELOCITY2>();
    auto u_pZp = cellADE.neighbor({0,0,1}).template getField<descriptors::VELOCITY2>();
    auto u_pZm = cellADE.neighbor({0,0,-1}).template getField<descriptors::VELOCITY2>();
    auto u_f = cellADE.template getField<descriptors::VELOCITY>();
    T por = (1. - cellADE.computeRho());
    T rho { };
    cellNSE.computeRhoU(rho, u_f.data());
 
    T relU[3] = {0.};
    relU[0] = (u_f[0] - u_p[0]) * convVel;
    relU[1] = (u_f[1] - u_p[1]) * convVel;
    relU[2] = (u_f[2] - u_p[2]) * convVel;
    T norm = util::sqrt(relU[0]*relU[0] + relU[1]*relU[1] + relU[2]*relU[2]);
    T Re_p = norm * particleDiam / visc;
    T C_D = 0.;
    if(Re_p != 0.){
    if(Re_p <= 1000.){
      C_D = 24.*(1. + 0.15*util::pow(Re_p,0.687))/Re_p;
    }else{
      C_D = 0.44;
    }}
    T F_D[3] = {0.};
    //F_D[0] = C_D * 3./4. * rho * convDens * norm * relU[0] / particleDiam / partDens;
    //F_D[1] = C_D * 3./4. * rho * convDens * norm * relU[1] / particleDiam / partDens;
    //F_D[2] = C_D * 3./4. * rho * convDens * norm * relU[2] / particleDiam / partDens;
    //dragCoeff = (9.*converter_.getPhysViscosity()*converter_.getPhysDensity()*converter_.getConversionFactorTime()) / (2.*pRho_*pRadius_*pRadius_);
    F_D[0] = 9.*visc*convDens/partDens/particleDiam*relU[0];
    F_D[1] = 9.*visc*convDens/partDens/particleDiam*relU[1];
    F_D[2] = 9.*visc*convDens/partDens/particleDiam*relU[2];
 
    T mass = partDens * 3.14/4. * particleDiam * particleDiam;
    if(cellADE.computeRho() > 0.001){
    u_p[0] += dt * (F_D[0]/*mass*/ + earthAcc[0]*(partDens - convDens)/partDens) / convVel - u_p * 0.5 * (u_pXp - u_pXm);
    u_p[1] += dt * (F_D[1]/*mass*/ + earthAcc[1]*(partDens - convDens)/partDens) / convVel - u_p * 0.5 * (u_pYp - u_pYm);
    u_p[2] += dt * (F_D[2]/*mass*/ + earthAcc[2]*(partDens - convDens)/partDens) / convVel - u_p * 0.5 * (u_pZp - u_pZm);
    cellADE.template setField<descriptors::VELOCITY>(u_p);
    }
 
    T porPlusX = (1. - cellADE.neighbor({1,0,0}).computeRho());
    T porMinusX = (1. - cellADE.neighbor({-1,0,0}).computeRho());
    T porPlusY = (1. - cellADE.neighbor({0,1,0}).computeRho());
    T porMinusY = (1. - cellADE.neighbor({0,-1,0}).computeRho());
    T porPlusZ = (1. - cellADE.neighbor({0,0,1}).computeRho());
    T porMinusZ = (1. - cellADE.neighbor({0,0,-1}).computeRho());
 
    T porDX2 = porPlusX -2.*por + porMinusX;
    T porDY2 = porPlusY -2.*por + porMinusY;
    T porDZ2 = porPlusZ -2.*por + porMinusZ;
 
    int choice = 0;
    if( porDX2 != 0.) choice += 1;
    if( porDY2 != 0.) choice += 1;
    if( porDZ2 != 0.) choice += 1;
 
    T coeff = (choice == 1) *1./4. + (choice == 2) * 1./6. + (choice == 3) * 5./36.;
    T pressCorr = por + coeff * (porDX2 + porDY2 + porDZ2);
    if(pressCorr < 0.001){ pressCorr = 1.; }
    cellNSE.template setField<descriptors::PRESSCORR>(pressCorr);
    T press = rho / pressCorr / descriptors::invCs2<T,DESCRIPTOR>();
    T pressCorrForce[3] = {0.};
    pressCorrForce[0] = press * 0.5 * (porPlusX - porMinusX);
    pressCorrForce[1] = press * 0.5 * (porPlusY - porMinusY);
    pressCorrForce[2] = press * 0.5 * (porPlusZ - porMinusZ);
 
    T force[3] = {0.};
    if(cellADE.computeRho() > 0.001){
    force[0] = pressCorrForce[0] / rho - (1. - por) * por * F_D[0] / convForce * convMass * partDens / rho + por * earthAcc[0] / convVel * dt;
    force[1] = pressCorrForce[1] / rho - (1. - por) * por * F_D[1] / convForce * convMass * partDens / rho + por * earthAcc[1] / convVel * dt;
    force[2] = pressCorrForce[2] / rho - (1. - por) * por * F_D[2] / convForce * convMass * partDens / rho + por * earthAcc[2] / convVel * dt;
    }
    cellNSE.template setField<descriptors::FORCE>(force);
  }