LPTEllipsoidal< nDim > Class Template Reference

#include <lptellipsoidal.h>

Inheritance diagram for LPTEllipsoidal< nDim >:

Collaboration diagram for LPTEllipsoidal< nDim >:

Public Member Functions
	LPTEllipsoidal ()
	constructor of ellipsoidal particles More...
	~LPTEllipsoidal () override=default
	LPTEllipsoidal (const LPTEllipsoidal &)=default
LPTEllipsoidal &	operator= (const LPTEllipsoidal &)=default
void	motionEquation () override
void	energyEquation () override
void	coupling () override
	single particle coupling terms More...
void	advanceParticle () override
	advance to new timeStep More...
void	resetWeights () override
void	calculateMajorAxisOrientation (MFloat *majorAxis)
	Calculate the orientation of the major axis in the world coordinate system. More...
MFloat	dragFactor (const MFloat partRe, const MFloat beta=1, const MFloat inclinationAngle=0, const MFloat K0=0, const MFloat K2=0)
	Calculates the drag factor for current particle with given Re number, aspect ratio and inclination angle. More...
MFloat	liftFactor (const MFloat partRe, const MFloat beta, const MFloat inclinationAngle, const MFloat K0, const MFloat K2)
	Calculates the lift factor for current particle with given Re number, aspect ratio, and inclination angle. More...
MFloat	torqueFactor (const MFloat partRe, const MFloat beta, const MFloat inclinationAngle)
	Calculates the torque factor for current particle with given Re number, aspect ratio and inclination angle. More...
void	setParticleFluidVelocities (MInt cellId, MFloat *position, MFloat *quaternions, MFloat *fluidVel, MFloat *fluidVelGrad)
MFloat	maxParticleRadius () const
	Returns the largest radius. More...
MFloat	minParticleRadius () const
	Returns the smallest radius. More...
MFloat	particleRadius (MInt direction) const
	Returns the radius in direction d. More...
MFloat	equivalentRadius () const
	Returns the equivalent radius. More...
MFloat	equivalentDiameter () const
	Returns the equivalent diameter. More...
MFloat	particleVolume () const
	Calculates the current volume. More...
MFloat	particleMass () const
	Calculates the current mass. More...
MFloat	magRelVel (const MFloat *const velocity1, const MFloat *const velocity2) const
	Calculates the magnitude of the relative velocity of the two given velocity vectors. More...
MFloat	particleRe (const MFloat velocity, const MFloat density, const MFloat dynamicViscosity)
	Calculates the particle Reynolds number. More...
MFloat	fParticleRelTime (const MFloat dynamicViscosity)
	Calculates the reciprocal of the particle relaxation time. More...
MFloat	effectiveWallCollisionRadius () const override
	Returns the relevant particle radius for the wall collision. More...
void	initEllipsoialProperties ()
void	initVelocityAndDensity ()
Public Member Functions inherited from LPTBase< nDim >
virtual	~LPTBase ()=default
BitsetType::reference	hasProperty (const LptBaseProperty p)
	Accessor for properties. More...
MBool	hasProperty (const LptBaseProperty p) const
	Accessor for properties (const version). More...
BitsetType::reference	isWindow ()
MBool	isWindow () const
BitsetType::reference	reqSend ()
MBool	reqSend () const
BitsetType::reference	reqBroadcast ()
MBool	reqBroadcast () const
BitsetType::reference	isInvalid ()
MBool	isInvalid () const
BitsetType::reference	hasCollided ()
MBool	hasCollided () const
BitsetType::reference	hadWallColl ()
MBool	hadWallColl () const
BitsetType::reference	firstStep ()
MBool	firstStep () const
BitsetType::reference	toBeDeleted ()
MBool	toBeDeleted () const
BitsetType::reference	wasSend ()
MBool	wasSend () const
BitsetType::reference	toBeRespawn ()
MBool	toBeRespawn () const
BitsetType::reference	fullyEvaporated ()
MBool	fullyEvaporated () const
void	getNghbrList (std::vector< MInt > &, const MInt)
template<MInt a, MInt b>
void	interpolateAndCalcWeights (const MInt cellId, const MFloat *const x, MFloat *const result, std::vector< MFloat > &weight, MFloat *const gradientResult=nullptr)
virtual void	energyEquation ()=0
virtual void	coupling ()=0
virtual void	motionEquation ()=0
virtual void	advanceParticle ()=0
virtual void	resetWeights ()=0
virtual void	particleWallCollision ()
	particle-wall collision step More...
virtual void	wallParticleCollision ()
	wall-particle collision step More...
void	initProperties ()
void	checkCellChange (const MFloat *oldPosition, const MBool allowHaloNonLeaf=false)
	Checks whether position is within cell with number cellId if not, cellId is updated. More...
void	updateProperties (const MBool init=true)
	Update particle properties. More...
virtual MFloat	effectiveWallCollisionRadius () const

Public Attributes
std::array< MFloat, nDim >	m_angularVel {}
	particle angular velocity More...
std::array< MFloat, nDim >	m_angularAccel {}
	particle angular acceleration More...
std::array< MFloat, nDim >	m_fluidVel {}
	fluid velocity More...
MFloat	m_fluidDensity = std::numeric_limits<MFloat>::quiet_NaN()
	fluid density More...
std::array< MFloat, nDim >	m_oldAccel {}
	particle acceleration of the last time step More...
std::array< MFloat, nDim >	m_oldAngularVel {}
	particle angular velocity of the last time step More...
std::array< MFloat, nDim >	m_oldAngularAccel {}
	particle angular acceleration of the last time step More...
std::array< MFloat, nDim >	m_oldFluidVel {}
	fluid velocity of the last time step More...
std::array< MFloat, nDim *nDim >	m_velocityGradientFluid {}
	gradient of fluid velocity in the particle fixed coordinate system More...
MFloat	m_oldFluidDensity = std::numeric_limits<MFloat>::quiet_NaN()
	old fluid density More...
MFloat	m_aspectRatio = std::numeric_limits<MFloat>::quiet_NaN()
	aspect ratio a/c of ellipsoidal: 0 < oblate < 1, =1 spherical, >1 prolate More...
MFloat	m_semiMinorAxis = std::numeric_limits<MFloat>::quiet_NaN()
	semi minor axis a More...
std::array< MFloat, 4 >	m_shapeParams {std::numeric_limits<MFloat>::quiet_NaN()}
	Shape parameters for ellipsoid. More...
std::array< MFloat, 4 >	m_quaternion {}
std::array< MFloat, 4 >	m_oldQuaternion {}
MFloat	m_temperature = std::numeric_limits<MFloat>::quiet_NaN()
	temperature More...
MFloat	m_fluidVelMag = std::numeric_limits<MFloat>::quiet_NaN()
	fluid velocity magnitude More...
MFloat	m_heatFlux = std::numeric_limits<MFloat>::quiet_NaN()
	heat flux to cell More...
MInt	m_particleMajorAxis = 2
std::vector< MFloat >	m_redistWeight {}
Public Attributes inherited from LPTBase< nDim >
std::array< MFloat, nDim >	m_position {}
std::array< MFloat, nDim >	m_velocity {}
std::array< MFloat, nDim >	m_accel {}
std::array< MFloat, nDim >	m_oldPos {}
	particle position of the last time step More...
std::array< MFloat, nDim >	m_oldVel {}
	particle velocity of the last time step More...
MFloat	m_densityRatio = -2
MLong	m_partId = -2
MInt	m_cellId = -2
MInt	m_oldCellId = -2
MFloat	m_creationTime = std::numeric_limits<MFloat>::quiet_NaN()
	creation time modifier More...
BitsetType	m_properties
std::vector< MInt >	m_neighborList

Static Public Attributes
static constexpr MInt	s_floatElements = 7 + 8 * nDim + 4 * 2
static constexpr MInt	s_intElements = 0
static MFloat	s_lengthFactor = F1
static MFloat	s_Re = F1
static MFloat	s_Pr = F1
static MFloat	s_Sc = F1
static MFloat	s_We
static std::array< MFloat, nDim >	s_Frm {}
Static Public Attributes inherited from LPTBase< nDim >
static MInt	s_interpolationOrder = 0
static MInt	s_interpolationMethod = 0
static MFloat	s_distFactorImp = 0.0
static LPT< nDim > *	s_backPtr = nullptr
static constexpr MInt	s_floatElements = 3 * nDim + 1
static constexpr MInt	s_intElements = 5

Private Member Functions
void	initShapeParams ()
	Calcultes shape parameters for current particle (see Siewert, 2014) More...
void	initEqDiameter ()
void	interpolateVelocityAndDensity (const MInt cellId, const MFloat *const position, MFloat *const velocity, MFloat &density, std::vector< MFloat > &weights)
void	interpolateVelocityAndVelocitySlopes (const MInt cellId, const MFloat *const position, MFloat *const velocity, MFloat *const velocityGradient, std::vector< MFloat > &weights)
void	momentumCoupling ()
	add the momentum flux from the particle to the cell/all surrounding cells More...
void	heatCoupling ()
	add the heat flux from the particle to the cell/all surrounding cells More...
void	massCoupling ()
	add the mass flux from the particle to the cell/all surrounding cells More...

Private Attributes
MFloat	m_eqDiameter = std::numeric_limits<MFloat>::quiet_NaN()
	equivalent particle diameter More...

Friends
template<MInt >
std::ostream &	operator<< (std::ostream &os, const LPTEllipsoidal &m)

Additional Inherited Members
Public Types inherited from LPTBase< nDim >
using	BitsetType = maia::lpt::baseProperty::BitsetType

Detailed Description

template<MInt nDim>
class LPTEllipsoidal< nDim >

Definition at line 21 of file lptellipsoidal.h.

Constructor & Destructor Documentation

LPTEllipsoidal() [1/2]

template<MInt nDim>

LPTEllipsoidal< nDim >::LPTEllipsoidal

Author

Laurent Andre

Definition at line 40 of file lptellipsoidal.cpp.

                                     {
  const auto nan = std::numeric_limits<MFloat>::quiet_NaN();
  for(MInt i = 0; i < nDim; i++) {
    m_fluidVel[i] = nan;
    m_oldFluidVel[i] = nan;
    m_velocity[i] = nan;
    m_oldVel[i] = nan;
    m_position[i] = nan;
    m_oldPos[i] = nan;
    m_accel[i] = nan;
    m_oldAccel[i] = nan;
    m_angularVel[i] = nan;
    m_oldAngularVel[i] = nan;
    m_angularAccel[i] = nan;
    m_oldAngularAccel[i] = nan;
  }
  for(MInt i = 0; i < 4; i++) {
    m_quaternion[i] = nan;
    m_oldQuaternion[i] = nan;
  }
}

~LPTEllipsoidal()

template<MInt nDim>

LPTEllipsoidal< nDim >::~LPTEllipsoidal ( )

overridedefault

LPTEllipsoidal() [2/2]

template<MInt nDim>

LPTEllipsoidal< nDim >::LPTEllipsoidal ( const LPTEllipsoidal< nDim > &  )

default

Member Function Documentation

advanceParticle()

template<MInt nDim>

void LPTEllipsoidal< nDim >::advanceParticle

overridevirtual

Author

Tim Wegmann

Implements LPTBase< nDim >.

Definition at line 68 of file lptellipsoidal.cpp.

                                           {
  firstStep() = false;
  hasCollided() = false;
  hadWallColl() = false;
  for(MInt i = 0; i < nDim; i++) {
    ASSERT(!isnan(m_fluidVel[i]), "");
    m_oldFluidVel[i] = m_fluidVel[i];
  }
  ASSERT(m_fluidDensity > 0, "");
  m_oldFluidDensity = m_fluidDensity;
}

calculateMajorAxisOrientation()

template<MInt nDim>

void LPTEllipsoidal< nDim >::calculateMajorAxisOrientation

(

MFloat *

majorAxis

)

Parameters

majorAxis

- output array to store the major axis

Definition at line 896 of file lptellipsoidal.cpp.

                                                                          {
  std::array<MFloat, nDim> majorAxisParticleFixed{}; // major axis in particle fixed coordinate system
  for(MInt i = 0; i < nDim; i++) {
    i == m_particleMajorAxis ? majorAxisParticleFixed[i] = maxParticleRadius() : majorAxisParticleFixed[i] = F0;
  }
  MFloatScratchSpace R(3, 3, AT_, "R");
  maia::math::computeRotationMatrix(R, &m_quaternion[0]);
  maia::math::matrixVectorProductTranspose(majorAxis, R, majorAxisParticleFixed.begin());
}

coupling()

template<MInt nDim>

void LPTEllipsoidal< nDim >::coupling

overridevirtual

Author

Tim Wegmann

Implements LPTBase< nDim >.

Definition at line 85 of file lptellipsoidal.cpp.

                                    {
  if(m_cellId < 0) {
    mTerm(1, AT_, "coupling of invalid particle?");
  }
 
  if(!s_backPtr->m_heatCoupling && !s_backPtr->m_evaporation) {
    m_redistWeight.clear();
    interpolateVelocityAndDensity(m_cellId, m_position.data(), m_fluidVel.data(), m_fluidDensity, m_redistWeight);
  }
 
  // check that all variables have been exchanged/set correctly
  for(MInt i = 0; i < nDim; i++) {
    ASSERT(!isnan(m_oldFluidVel[i]), "");
  }
  ASSERT(m_oldFluidDensity > 0, "");
 
  // check all source terms
#ifdef LPT_DEBUG
  for(MInt i = 0; i < nDim; i++) {
    if(std::isnan(m_accel[i])) {
      cerr << "Invalid motion equation! " << m_partId << " " << m_position[0] << " " << m_position[1] << " "
           << m_position[2] << endl;
      isInvalid() = true;
      return;
    }
    if(std::isnan(m_velocity[i])) {
      cerr << "Invalid motion equation2! " << m_partId << " " << m_position[0] << " " << m_position[1] << " "
           << m_position[2] << " " << m_densityRatio << " " << hadWallColl() << endl;
      isInvalid() = true;
      return;
    }
  }
#endif
 
  if(s_backPtr->m_momentumCoupling) momentumCoupling();
  if(s_backPtr->m_heatCoupling) heatCoupling();
  if(s_backPtr->m_massCoupling) massCoupling();
}

dragFactor()

template<MInt nDim>

MFloat LPTEllipsoidal< nDim >::dragFactor	(	const MFloat	partRe,
		const MFloat	beta = 1,
		const MFloat	inclinationAngle = 0,
		const MFloat	K0 = 0,
		const MFloat	K2 = 0
	)

Author

Laurent Andre

Date

August 2022

Definition at line 741 of file lptellipsoidal.cpp.

                                                                          {
  MFloat result = NAN;
  switch(s_backPtr->m_dragModelType) {
    case 0: { // drag force switched off
      result = F0;
      break;
    }
    case 1: { // linear Stokes drag
      result = F1;
      break;
    }
    case 2: { // nonlinear drag according to Schiller & Naumann
      result = 1.0 + 0.15 * pow(partRe, 0.687);
      break;
    }
    case 3: { // new correlations by Konstantin (2020)
      std::array<MFloat, 8> dCorrs{-0.007, 1.0, 1.17, -0.07, 0.047, 1.14, 0.7, -0.008};
      MFloat fd0 = F1;
      MFloat fd90 = F1;
      fd0 = 1.0 + 0.15 * pow(partRe, 0.687)
            + dCorrs[0] * pow(std::log(beta), dCorrs[1]) * pow(partRe, dCorrs[2] + dCorrs[3] * std::log(beta));
      fd90 = 1.0 + 0.15 * pow(partRe, 0.687)
             + dCorrs[4] * pow(log(beta), dCorrs[5]) * pow(partRe, dCorrs[6] + dCorrs[7] * log(beta));
 
      MFloat fd0K2 = fd0 * K2;
      MFloat f90K0 = fd90 * K0;
      result = fd0K2 + (f90K0 - fd0K2) * POW2(sin(inclinationAngle));
      break;
    }
    default: {
      mTerm(1, AT_, "Unknown particle drag method");
    }
  }
  return result;
}

effectiveWallCollisionRadius()

template<MInt nDim>

MFloat LPTEllipsoidal< nDim >::effectiveWallCollisionRadius ( ) const

inlineoverridevirtual

Reimplemented from LPTBase< nDim >.

Definition at line 186 of file lptellipsoidal.h.

{ return equivalentRadius(); }

energyEquation()

template<MInt nDim>

void LPTEllipsoidal< nDim >::energyEquation

overridevirtual

Implements LPTBase< nDim >.

Definition at line 885 of file lptellipsoidal.cpp.

                                          {
  mTerm(1, AT_, "Not yet implemented");
}

equivalentDiameter()

template<MInt nDim>

MFloat LPTEllipsoidal< nDim >::equivalentDiameter ( ) const

inline

Definition at line 156 of file lptellipsoidal.h.

{ return m_eqDiameter; }

equivalentRadius()

template<MInt nDim>

MFloat LPTEllipsoidal< nDim >::equivalentRadius ( ) const

inline

Definition at line 153 of file lptellipsoidal.h.

{ return F1B2 * m_eqDiameter; }

fParticleRelTime()

template<MInt nDim>

MFloat LPTEllipsoidal< nDim >::fParticleRelTime ( const MFloat  dynamicViscosity )

inline

Definition at line 181 of file lptellipsoidal.h.

                                                                {
    return 18.0 * dynamicViscosity / (m_eqDiameter * m_eqDiameter * s_backPtr->m_material->density());
  }

heatCoupling()

template<MInt nDim>

void LPTEllipsoidal< nDim >::heatCoupling

private

Author

Sven Berger, Tim Wegmann

Definition at line 179 of file lptellipsoidal.cpp.

                                        {
  mTerm(1, AT_, "Heat coupling not tested for ellipsoidal particles");
}

initEllipsoialProperties()

template<MInt nDim>

void LPTEllipsoidal< nDim >::initEllipsoialProperties ( )

inline

Definition at line 188 of file lptellipsoidal.h.

                                         {
    ASSERT(!std::isnan(m_aspectRatio), "Aspect Ratio not set!");
    ASSERT(!std::isnan(m_semiMinorAxis), "Semi Minor Axis not set!");
    initEqDiameter();
    initShapeParams();
  }

initEqDiameter()

template<MInt nDim>

void LPTEllipsoidal< nDim >::initEqDiameter ( )

inlineprivate

Definition at line 212 of file lptellipsoidal.h.

{ m_eqDiameter = 2 * m_semiMinorAxis * pow(m_aspectRatio, F1B3); }

initShapeParams()

template<MInt nDim>

void LPTEllipsoidal< nDim >::initShapeParams

private

Author

Laurent Andre

Date

August 2022

Definition at line 854 of file lptellipsoidal.cpp.

                                           {
  MFloat beta = m_aspectRatio;
  MFloat beta2 = beta * beta;
  if(fabs(beta - F1) < 1e-14) {
    // spherical
    m_shapeParams[0] = F2 * POW2(m_semiMinorAxis);
    m_shapeParams[1] = F2B3;
    m_shapeParams[2] = F2B3;
    m_shapeParams[3] = F2B3;
  } else if(beta > F1) {
    // prolate
    MFloat kappa = log((beta - sqrt(beta2 - F1)) / (beta + sqrt(beta2 - F1)));
    m_shapeParams[0] = -beta * POW2(m_semiMinorAxis) * kappa / sqrt(beta2 - F1);
    m_shapeParams[1] = beta2 / (beta2 - F1) + beta * kappa / (F2 * pow(beta2 - F1, F3B2));
    m_shapeParams[2] = beta2 / (beta2 - F1) + beta * kappa / (F2 * pow(beta2 - F1, F3B2));
    m_shapeParams[3] = -F2 / (beta2 - F1) - beta * kappa / (pow(beta2 - F1, F3B2));
  } else if(beta < F1) {
    // oblate
    MFloat kappa = F2 * atan2(beta, sqrt(F1 - beta2));
    m_shapeParams[0] = beta * POW2(m_semiMinorAxis) * (PI - kappa) / sqrt(F1 - beta2);
    m_shapeParams[1] = -beta * (kappa - PI + F2 * beta * sqrt(F1 - beta2)) / (F2 * pow(F1 - beta2, F3B2));
    m_shapeParams[1] = -beta * (kappa - PI + F2 * beta * sqrt(F1 - beta2)) / (F2 * pow(F1 - beta2, F3B2));
    m_shapeParams[1] = (beta * kappa - beta * PI + F2 * sqrt(F1 - beta2)) / (pow(F1 - beta2, F3B2));
  } else {
    // general triaxial ellipsoid
    mTerm(1, "Generalize here");
  }
}

initVelocityAndDensity()

template<MInt nDim>

void LPTEllipsoidal< nDim >::initVelocityAndDensity ( )

inline

Definition at line 195 of file lptellipsoidal.h.

                                {
    std::array<MFloat, nDim + 1> result{};
    std::vector<MFloat> weights;
    this->template interpolateAndCalcWeights<0, nDim + 1>(m_cellId, m_position.data(), result.data(), weights);
    for(MInt i = 0; i < nDim; i++) {
      m_oldFluidVel[i] = result[i];
      m_fluidVel[i] = result[i];
    }
    m_oldFluidDensity = result[nDim];
    m_fluidDensity = result[nDim];
  }

interpolateVelocityAndDensity()

template<MInt nDim>

void LPTEllipsoidal< nDim >::interpolateVelocityAndDensity ( const MInt  cellId, const MFloat *const  position, MFloat *const  velocity, MFloat &  density, std::vector< MFloat > &  weights  )

inlineprivate

Definition at line 214 of file lptellipsoidal.h.

                                                                                  {
    std::array<MFloat, nDim + 1> result{};
    this->template interpolateAndCalcWeights<0, nDim + 1>(cellId, position, result.data(), weights);
    for(MInt n = 0; n < nDim; n++) {
      velocity[n] = result[n];
    }
    density = result[nDim];
  }

interpolateVelocityAndVelocitySlopes()

template<MInt nDim>

void LPTEllipsoidal< nDim >::interpolateVelocityAndVelocitySlopes ( const MInt  cellId, const MFloat *const  position, MFloat *const  velocity, MFloat *const  velocityGradient, std::vector< MFloat > &  weights  )

inlineprivate

Definition at line 224 of file lptellipsoidal.h.

                                                                                                        {
    std::array<MFloat, nDim> result{};
    std::array<MFloat, nDim * nDim> gradientResult{};
    this->template interpolateAndCalcWeights<0, nDim>(cellId, position, result.data(), weights, gradientResult.data());
    for(MInt n = 0; n < nDim; n++) {
      velocity[n] = result[n];
    }
    for(MInt n = 0; n < nDim * nDim; n++) {
      velocityGradient[n] = gradientResult[n];
    }
  }

liftFactor()

template<MInt nDim>

MFloat LPTEllipsoidal< nDim >::liftFactor	(	const MFloat	partRe,
		const MFloat	beta,
		const MFloat	inclinationAngle,
		const MFloat	K0,
		const MFloat	K2
	)

Author

Laurent Andre

Date

August 2022

Definition at line 784 of file lptellipsoidal.cpp.

                                                                          {
  MFloat result = NAN;
  switch(s_backPtr->m_liftModelType) {
    case 0: { // lift force switched off
      result = F0;
      break;
    }
    case 1: { // linear
      result = F1;
      break;
    }
    case 2: { // new correlations by Konstantin (2020)
      std::array<MFloat, 6> lCorrs{0.01, 0.86, 1.77, 0.34, 0.88, -0.05};
      MFloat flMax = F1;
      MFloat flShift = F1;
      if(partRe > 1) flShift += lCorrs[0] * pow(log(beta), lCorrs[1]) * pow(log(partRe), lCorrs[2]);
      MFloat psi = F1B2 * M_PI * pow(inclinationAngle / M_PI * F2, flShift);
      flMax = F1 + lCorrs[3] * pow(partRe, lCorrs[4] + lCorrs[5]);
      MFloat ClMax = F1B2 * flMax * (K0 - K2);
      MFloat Cl = F2 * sin(psi) * cos(psi) * ClMax;
      result = Cl;
      break;
    }
    default: {
      mTerm(1, AT_, "Unknown particle drag method");
    }
  }
  return result;
}

magRelVel()

template<MInt nDim>

MFloat LPTEllipsoidal< nDim >::magRelVel ( const MFloat *const  velocity1, const MFloat *const  velocity2  ) const

inline

Definition at line 167 of file lptellipsoidal.h.

                                                                                              {
    MFloat magnitude = 0.0;
    for(MInt i = 0; i < nDim; i++) {
      magnitude += POW2(velocity1[i] - velocity2[i]);
    }
    return sqrt(magnitude);
  }

massCoupling()

template<MInt nDim>

void LPTEllipsoidal< nDim >::massCoupling

private

Author

Sven Berger, Tim Wegmann

Definition at line 129 of file lptellipsoidal.cpp.

                                        {
  mTerm(1, AT_, "Mass coupling not implemented for ellipsoidal particles");
}

maxParticleRadius()

template<MInt nDim>

MFloat LPTEllipsoidal< nDim >::maxParticleRadius ( ) const

inline

Definition at line 141 of file lptellipsoidal.h.

{ return m_semiMinorAxis * m_aspectRatio; }

minParticleRadius()

template<MInt nDim>

MFloat LPTEllipsoidal< nDim >::minParticleRadius ( ) const

inline

Definition at line 143 of file lptellipsoidal.h.

{ return m_semiMinorAxis; }

momentumCoupling()

template<MInt nDim>

void LPTEllipsoidal< nDim >::momentumCoupling

private

Author

Sven Berger, Tim Wegmann

Definition at line 138 of file lptellipsoidal.cpp.

                                            {
  const MFloat relT = (s_backPtr->m_timeStep - m_creationTime) / s_backPtr->m_timeStep;
 
  const MFloat mass = particleMass();
 
  array<MFloat, nDim> force{};
 
  for(MInt i = 0; i < nDim; i++) {
    // momentum flux due to drag force
    const MFloat invDensityRatio = (abs(m_densityRatio) <= MFloatEps) ? 1.0 : 1.0 / m_densityRatio;
    force[i] = mass * (m_accel[i] - ((1.0 - invDensityRatio) * s_Frm[i])) * relT;
  }
 
  // work done by the momentum
  const MFloat work = std::inner_product(force.begin(), force.end(), m_oldVel.begin(), 0.0);
 
 
  if(!s_backPtr->m_couplingRedist) {
    for(MInt i = 0; i < nDim; i++) {
      s_backPtr->a_momentumFlux(m_cellId, i) += (force[i]);
    }
    s_backPtr->a_workFlux(m_cellId) += work;
 
  } else {
    ASSERT(m_redistWeight.size() == m_neighborList.size(),
           "Invalid " + to_string(m_redistWeight.size()) + "!=" + to_string(m_neighborList.size()));
    for(MInt n = 0; n < (signed)m_neighborList.size(); n++) {
      // a weightHeat is used to distribute the source term to multiple cells...
      for(MInt i = 0; i < nDim; i++) {
        s_backPtr->a_momentumFlux(m_neighborList[n], i) += m_redistWeight.at(n) * (force[i]);
      }
      s_backPtr->a_workFlux(m_neighborList[n]) += m_redistWeight.at(n) * work;
    }
  }
}

motionEquation()

template<MInt nDim>

void LPTEllipsoidal< nDim >::motionEquation

overridevirtual

1. Predictor step
2. get velocity, density and velocity magnitude at the predicted position
3. corrector step

Implements LPTBase< nDim >.

Definition at line 185 of file lptellipsoidal.cpp.

                                          {
  const MFloat dt = s_backPtr->m_timeStep - m_creationTime;
  ASSERT(dt > 0, "Timestep < 0");
 
  // NOTE: at this point old means 1 time-Step ago and the current values will be updated below!
  m_oldPos = m_position;
  m_oldVel = m_velocity;
  m_oldAccel = m_accel;
  m_oldAngularVel = m_angularVel;
  m_oldAngularAccel = m_angularAccel;
  m_oldQuaternion = m_quaternion;
  m_oldCellId = m_cellId;
 
  if(firstStep()) {
    for(MInt i = 0; i < nDim; i++) {
      m_accel[i] = s_Frm[i];
      m_oldAccel[i] = s_Frm[i];
    }
  }
 
  for(MInt i = 0; i < nDim; i++) {
    ASSERT(!isnan(m_oldFluidVel[i]), "");
  }
  ASSERT(m_oldFluidDensity > 0, "");
 
  // reduce variance between runs (round to 9 decimal places)
  // TODO-timw labels:LPT,totest check if rounding is still necessary?
  const MFloat invDensityRatio =
      F1 / (static_cast<MFloat>(ceil(s_backPtr->m_material->density() / m_oldFluidDensity * 1E9)) / 1E9);
 
  m_densityRatio = 1 / invDensityRatio;
 
  array<MFloat, nDim> predictedPos{};
  array<MFloat, nDim> predictedVel{};
  array<MFloat, nDim> predictedAngularVel{};
  array<MFloat, nDim + 1> predictedQuaternion{};
 
  MFloat q[3];
  MFloat q0[3];
  MFloat tq[3];
  MFloat tq0[3];
  MFloat w = m_oldQuaternion[0];
  MFloat x = m_oldQuaternion[1];
  MFloat y = m_oldQuaternion[2];
  MFloat z = m_oldQuaternion[3];
  for(MInt i = 0; i < 3; i++) {
    q0[i] = m_oldAngularVel[i];
  }
 
  // predict position, velocity, orientation and angular velocity
  for(MInt i = 0; i < nDim; i++) {
    predictedVel[i] = m_oldVel[i] + dt * m_oldAccel[i];
    predictedPos[i] = m_oldPos[i] + dt * (m_oldVel[i]) + F1B2 * POW2(dt) * (m_oldAccel[i]);
  }
  predictedQuaternion[0] = w + F1B2 * dt * (-x * q0[0] - y * q0[1] - z * q0[2]);
  predictedQuaternion[1] = x + F1B2 * dt * (w * q0[0] - z * q0[1] + y * q0[2]);
  predictedQuaternion[2] = y + F1B2 * dt * (z * q0[0] + w * q0[1] - x * q0[2]);
  predictedQuaternion[3] = z + F1B2 * dt * (-y * q0[0] + x * q0[1] + w * q0[2]);
  for(MInt i = 0; i < 3; i++) {
    predictedAngularVel[i] = m_oldAngularVel[i] + dt * m_oldAngularAccel[i];
  }
 
  //  find the predicted CellId at the predicted position around the previous cellId
  const MInt predictedCellId = s_backPtr->grid().findContainingLeafCell(&predictedPos[0], m_cellId);
 
  ASSERT(m_cellId > -1, "Invalid cell!");
  ASSERT(predictedCellId > -1, "Invalid cell!");
  ASSERT(s_backPtr->c_isLeafCell(predictedCellId), "Invalid cell!");
 
  vector<MFloat> predictedWeights(0);
  array<MFloat, nDim> fluidVel{};
  array<MFloat, nDim * nDim> fluidVelGradient{};
  array<MFloat, nDim * nDim> fluidVelGradientHat{};
  std::fill(fluidVel.begin(), fluidVel.end(), std::numeric_limits<MFloat>::quiet_NaN());
  std::fill(fluidVelGradient.begin(), fluidVelGradient.end(), std::numeric_limits<MFloat>::quiet_NaN());
  std::fill(fluidVelGradientHat.begin(), fluidVelGradientHat.end(), std::numeric_limits<MFloat>::quiet_NaN());
  interpolateVelocityAndVelocitySlopes(predictedCellId, predictedPos.data(), fluidVel.data(), fluidVelGradient.data(),
                                       predictedWeights);
 
  // calculate gradients according to principle axes of ellipsoid
  MFloatScratchSpace R(3, 3, AT_, "R");
  maia::math::computeRotationMatrix(R, &(m_quaternion[0]));
  for(MInt i = 0; i < nDim; i++) {
    maia::math::matrixVectorProduct(&fluidVelGradientHat[nDim * i], R, &fluidVelGradient[nDim * i]);
  }
 
  m_fluidVelMag = sqrt(std::inner_product(fluidVel.begin(), fluidVel.end(), fluidVel.begin(), F0));
 
  const MFloat T = s_backPtr->a_fluidTemperature(m_cellId);
 
  MFloat fluidViscosity = s_backPtr->m_material->dynViscosityFun(T);
 
  if(s_backPtr->m_dragModelType == 0) {
    // no drag or neglible small (are going to be deleted...)
    for(MInt i = 0; i < nDim; i++) {
      m_velocity[i] = predictedVel[i];
      m_position[i] = predictedPos[i];
      m_angularVel[i] = predictedAngularVel[i];
      m_quaternion[i] = predictedQuaternion[i];
      m_accel[i] = m_oldAccel[i];
      m_angularAccel[i] = m_oldAngularAccel[i];
    }
  } else {
    switch(s_backPtr->m_motionEquationType) {
      case 0: {
        // particle moves with fluid velocity, no rotational motion
        for(MInt i = 0; i < nDim; i++) {
          m_position[i] = m_oldPos[i] + F1B2 * dt * (fluidVel[i] + m_oldVel[i]) * s_lengthFactor;
          m_velocity[i] = fluidVel[i];
          m_accel[i] = (m_velocity[i] - m_oldVel[i]) / dt;
          m_angularVel[i] = F0;
          m_angularAccel[i] = F0;
        }
        break;
      }
      case 1: {
        // spherical case copied from lptspherical
        // serves for debugging to have a version independant of the strain, vorticity
        const MFloat relVel = magRelVel(&fluidVel[0], &predictedVel[0]);
 
        const MFloat Rep = particleRe(relVel, m_oldFluidDensity, fluidViscosity) * s_Re;
 
        const MFloat CD = dragFactor(Rep) / s_Re * fParticleRelTime(fluidViscosity);
 
        for(MInt i = 0; i < nDim; i++) {
          m_accel[i] = CD * (fluidVel[i] - predictedVel[i]) + (1.0 - invDensityRatio) * s_Frm[i];
          m_velocity[i] = m_oldVel[i] + dt * F1B2 * (m_accel[i] + m_oldAccel[i]);
          m_position[i] = m_oldPos[i] + dt * F1B2 * (m_velocity[i] + m_oldVel[i]) * s_lengthFactor;
        }
        break;
      }
      case 2: {
        // case derived for creeping flow
        const MFloat fdtB2 = F1B2 * dt;
        MFloatScratchSpace K(3, 3, AT_, "K");
        MFloatScratchSpace W(4, 4, AT_, "W");
        MFloatScratchSpace iW(4, 4, AT_, "iW");
        MFloatScratchSpace rhs(4, AT_, "rhs");
 
        MFloat vrel[3];
        MFloat vrelhat[3];
        MFloat tmp[3];
        MFloat vort[3];
        MFloat strain[3];
 
        const MFloat beta = m_aspectRatio;
        const MFloat beta2 = POW2(beta);
        const MFloat fTauP = fParticleRelTime(fluidViscosity) / s_Re;
 
        const MFloat linAccFac = (8.0 / 3.0) * pow(beta, F2B3) * fTauP;
        const MFloat angAccFac = (40.0 / 9.0) * pow(beta, F2B3) * fTauP;
 
        // integrate angular velocity
        const MInt maxit = 100;
        MFloat delta = F1;
        MInt it = 0;
 
        // set temporary variables for angular velocity (q) and acceleration (tq)
        for(MInt i = 0; i < 3; i++) {
          tq0[i] = m_oldAngularAccel[i];
          q0[i] = m_oldAngularVel[i];
          tq[i] = tq0[i];
          q[i] = q0[i];
        }
 
        // compute vorticity and strain
        for(MInt i = 0; i < 3; i++) {
          MInt id0 = (i + 1) % 3;
          MInt id1 = (id0 + 1) % 3;
          strain[i] = F1B2 * (fluidVelGradientHat[3 * id1 + id0] + fluidVelGradientHat[3 * id0 + id1]);
          vort[i] = F1B2 * (fluidVelGradientHat[3 * id1 + id0] - fluidVelGradientHat[3 * id0 + id1]);
        }
 
        W(0, 0) = F1 + fdtB2 * angAccFac / (m_shapeParams[1] + beta2 * m_shapeParams[3]);
        W(1, 1) = F1 + fdtB2 * angAccFac / (m_shapeParams[1] + beta2 * m_shapeParams[3]);
        W(2, 2) = F1 + fdtB2 * angAccFac / (F2 * m_shapeParams[1]);
        W(2, 0) = F0;
        W(2, 1) = F0;
 
        // Newton iterations
        while(delta > 1e-10 && it < maxit) {
          W(0, 1) = -fdtB2 * q[2] * (beta2 - F1) / (beta2 + F1);
          W(0, 2) = -fdtB2 * q[1] * (beta2 - F1) / (beta2 + F1);
          W(1, 0) = -fdtB2 * q[2] * (F1 - beta2) / (beta2 + F1);
          W(1, 2) = -fdtB2 * q[0] * (F1 - beta2) / (beta2 + F1);
 
          maia::math::invert(W, iW, 3, 3);
 
          // compute angular acceleration tq
          tq[0] = q[1] * q[2] * (beta2 - F1) / (beta2 + F1)
                  + angAccFac * (((F1 - beta2) / (F1 + beta2)) * strain[0] + vort[0] - q[0])
                        / (m_shapeParams[1] + beta2 * m_shapeParams[3]);
          tq[1] = q[0] * q[2] * (F1 - beta2) / (beta2 + F1)
                  + angAccFac * (((beta2 - F1) / (F1 + beta2)) * strain[1] + vort[1] - q[1])
                        / (m_shapeParams[1] + beta2 * m_shapeParams[3]);
          tq[2] = angAccFac * (vort[2] - q[2]) / (F2 * m_shapeParams[1]);
 
          for(MInt i = 0; i < 3; i++)
            rhs[i] = q[i] - q0[i] - fdtB2 * (tq0[i] + tq[i]);
 
          delta = F0;
          for(MInt i = 0; i < 3; i++) {
            const MFloat qq = q[i];
            for(MInt j = 0; j < 3; j++) {
              q[i] -= iW(i, j) * rhs(j);
            }
            delta = mMax(delta, fabs(q[i] - qq));
          }
          it++;
        }
 
        if(it >= maxit || it <= 0) {
          cerr << "Newton iterations did not converge " << m_partId << endl;
        }
 
        for(MInt i = 0; i < 3; i++) {
          m_angularVel[i] = q[i];
          m_angularAccel[i] = tq[i];
        }
 
        // integrate quaternions
        w = m_oldQuaternion[0];
        x = m_oldQuaternion[1];
        y = m_oldQuaternion[2];
        z = m_oldQuaternion[3];
 
        W(0, 0) = F0;
        W(0, 1) = -q[0];
        W(0, 2) = -q[1];
        W(0, 3) = -q[2];
        W(1, 0) = q[0];
        W(1, 1) = F0;
        W(1, 2) = q[2];
        W(1, 3) = -q[1];
        W(2, 0) = q[1];
        W(2, 1) = -q[2];
        W(2, 2) = F0;
        W(2, 3) = q[0];
        W(3, 0) = q[2];
        W(3, 1) = q[1];
        W(3, 2) = -q[0];
        W(3, 3) = F0;
 
        for(MInt i = 0; i < 4; i++) {
          for(MInt j = 0; j < 4; j++) {
            W(i, j) = -F1B2 * fdtB2 * W(i, j);
          }
          W(i, i) = F1;
        }
 
        rhs(0) = w + F1B2 * fdtB2 * (-x * q0[0] - y * q0[1] - z * q0[2]);
        rhs(1) = x + F1B2 * fdtB2 * (w * q0[0] - z * q0[1] + y * q0[2]);
        rhs(2) = y + F1B2 * fdtB2 * (z * q0[0] + w * q0[1] - x * q0[2]);
        rhs(3) = z + F1B2 * fdtB2 * (-y * q0[0] + x * q0[1] + w * q0[2]);
 
        maia::math::invert(W, iW, 4, 4);
 
        for(MInt i = 0; i < 4; i++) {
          m_quaternion[i] = F0;
          for(MInt j = 0; j < 4; j++) {
            m_quaternion[i] += iW(i, j) * rhs(j);
          }
        }
 
        MFloat abs = F0;
        for(MInt i = 0; i < 4; i++)
          abs += POW2(m_quaternion[i]);
        for(MInt i = 0; i < 4; i++)
          m_quaternion[i] /= sqrt(abs);
 
        // integrate linear motion
        K.fill(F0);
        K(0, 0) = F1 / (m_shapeParams[0] / POW2(m_semiMinorAxis) + m_shapeParams[1]);
        K(1, 1) = K(0, 0);
        K(2, 2) = F1 / (m_shapeParams[0] / POW2(m_semiMinorAxis) + beta2 * m_shapeParams[3]);
        maia::math::computeRotationMatrix(R, &(m_quaternion[0]));
 
        for(MInt i = 0; i < nDim; i++) {
          vrel[i] = fluidVel[i] - predictedVel[i];
        }
 
        maia::math::matrixVectorProduct(vrelhat, R, vrel); // principal axes frame relative velocity
        maia::math::matrixVectorProduct(tmp, K, vrelhat);
        maia::math::matrixVectorProductTranspose(vrel, R, tmp);
 
        for(MInt i = 0; i < nDim; i++) {
          m_accel[i] = linAccFac * vrel[i] + (1.0 - invDensityRatio) * s_Frm[i];
          m_velocity[i] = m_oldVel[i] + dt * F1B2 * (m_accel[i] + m_oldAccel[i]);
          m_position[i] = m_oldPos[i] + dt * F1B2 * (m_velocity[i] + m_oldVel[i]) * s_lengthFactor
                          + F1B4 * POW2(dt) * (m_accel[i] + m_oldAccel[i]);
        }
        break;
      }
      case 3: {
        // case using new correlations functions by Konstantin (2020)
        const MFloat fdtB2 = F1B2 * dt;
        MFloatScratchSpace K(3, 3, AT_, "K");
        MFloatScratchSpace W(4, 4, AT_, "W");
        MFloatScratchSpace iW(4, 4, AT_, "iW");
        MFloatScratchSpace rhs(4, AT_, "rhs");
 
        MFloat vrel[3];
        MFloat vrelhat[3];
        MFloat tmp[3];
        MFloat vort[3];
        MFloat strain[3];
 
        const MFloat beta = m_aspectRatio;
        const MFloat beta2 = POW2(beta);
        const MFloat eqRadius = F1B2 * m_eqDiameter;
        const MFloat fTauP = fParticleRelTime(fluidViscosity) / s_Re;
        const MFloat magRelVeloc = magRelVel(&fluidVel[0], &predictedVel[0]);
        const MFloat Rep = particleRe(magRelVeloc, m_oldFluidDensity, fluidViscosity) * s_Re;
 
        const MFloat angAccFac = (40.0 / 9.0) * pow(beta, F2B3) * fTauP;
        const MFloat pitchingTorqueFac = (15.0 / 8.0) * invDensityRatio * pow(beta, F2B3) / POW2(eqRadius);
        const MFloat addFacCd =
            48.0 * fluidViscosity / (s_Re * s_backPtr->m_material->density() * POW2(m_eqDiameter)) * pow(beta, F2B3);
        const MFloat momI[3] = {1.0 + POW2(beta), 1.0 + POW2(beta), 2.0}; // factors for moments of inertia
 
        maia::math::computeRotationMatrix(R, &(m_quaternion[0]));
 
        for(MInt i = 0; i < nDim; i++) {
          vrel[i] = fluidVel[i] - predictedVel[i];
        }
 
        maia::math::matrixVectorProduct(vrelhat, R, vrel); // principal axes frame relative velocity
        maia::math::matrixVectorProduct(tmp, K, vrelhat);
        maia::math::matrixVectorProductTranspose(vrel, R, tmp);
 
        const MFloat magRelVelHat = sqrt(POW2(vrelhat[0]) + POW2(vrelhat[1]) + POW2(vrelhat[2]));
 
        const auto nan = std::numeric_limits<MFloat>::quiet_NaN();
        MFloat inclinationAngle = nan;
        std::array<MFloat, 3> accHat = {nan, nan, nan};
        std::array<MFloat, 3> accDragHat = {nan, nan, nan};
        std::array<MFloat, 3> accLiftHat = {nan, nan, nan};
        std::array<MFloat, 3> dirDHat = {nan, nan, nan};
        std::array<MFloat, 3> dirTHat = {nan, nan, nan};
        std::array<MFloat, 3> dirLHat = {nan, nan, nan};
        std::array<MFloat, 3> dirDCrossZ = {nan, nan, nan};
        std::array<MFloat, 3> dirZHat{nan, nan, nan};
 
        // dirZHat
        for(MInt i = 0; i < nDim; i++) {
          i == m_particleMajorAxis ? dirZHat[i] = F1 : dirZHat[i] = F0;
        }
        // dirDhat
        for(MInt i = 0; i < nDim; i++) {
          dirDHat[i] = vrelhat[i] / magRelVelHat;
        }
        // dirDcrossZ = dirDhat x dirZhat
        maia::math::cross(&dirDHat[0], &dirZHat[0], &dirDCrossZ[0]);
        const MFloat dotProduct = std::inner_product(dirDHat.begin(), dirDHat.end(), dirZHat.begin(), F0);
        inclinationAngle = acos(fabs(dotProduct));
        MInt sgn = (dotProduct > 0 ? 1 : -1);
        for(MInt i = 0; i < nDim; i++) {
          dirTHat[i] = dirDCrossZ[i] / sqrt(POW2(dirDCrossZ[0]) + POW2(dirDCrossZ[1]) + POW2(dirDCrossZ[2])) * sgn;
        }
        // dirLHat = dirDhat x dirThat
        maia::math::cross(&dirDHat[0], &dirTHat[0], &dirLHat[0]);
 
        // integrate angular velocity
        const MInt maxit = 100;
        MFloat delta = F1;
        MInt it = 0;
 
        W(0, 0) = F1 + fdtB2 * angAccFac / (m_shapeParams[1] + beta2 * m_shapeParams[3]);
        W(1, 1) = F1 + fdtB2 * angAccFac / (m_shapeParams[1] + beta2 * m_shapeParams[3]);
        W(2, 2) = F1 + fdtB2 * angAccFac / (F2 * m_shapeParams[1]);
        W(2, 0) = F0;
        W(2, 1) = F0;
 
        // set temporary variables for angular velocity (q) and acceleration (tq)
        for(MInt i = 0; i < 3; i++) {
          tq0[i] = m_oldAngularAccel[i];
          q0[i] = m_oldAngularVel[i];
          tq[i] = tq0[i];
          q[i] = q0[i];
        }
 
        // compute pitching torque unless torque modelling is turned off
        MFloat pitch[3] = {F0};
        if(s_backPtr->m_torqueModelType > 0 && Rep > 1e-8) {
          MFloat CT = torqueFactor(Rep, beta, inclinationAngle);
          for(MInt i = 0; i < 3; i++) {
            pitch[i] = pitchingTorqueFac * POW2(magRelVelHat) * CT / momI[i] * dirTHat[i];
            if(pitchingTorqueFac <= F0 || CT <= F0 || momI[i] <= F0)
              cout << "Calc pitching torque " << pitchingTorqueFac << " " << CT << " " << momI[i] << endl;
          }
        }
 
        // compute vorticity and strain
        for(MInt i = 0; i < 3; i++) {
          MInt id0 = (i + 1) % 3;
          MInt id1 = (id0 + 1) % 3;
          strain[i] = F1B2 * (fluidVelGradientHat[3 * id1 + id0] + fluidVelGradientHat[3 * id0 + id1]);
          vort[i] = F1B2 * (fluidVelGradientHat[3 * id1 + id0] - fluidVelGradientHat[3 * id0 + id1]);
        }
 
        // Newton iterations
        while(delta > 1e-10 && it < maxit) {
          W(0, 1) = -fdtB2 * q[2] * (beta2 - F1) / (beta2 + F1);
          W(0, 2) = -fdtB2 * q[1] * (beta2 - F1) / (beta2 + F1);
          W(1, 0) = -fdtB2 * q[2] * (F1 - beta2) / (beta2 + F1);
          W(1, 2) = -fdtB2 * q[0] * (F1 - beta2) / (beta2 + F1);
 
          maia::math::invert(W, iW, 3, 3);
 
          tq[0] = q[1] * q[2] * (beta2 - F1) / (beta2 + F1)
                  + angAccFac * (((F1 - beta2) / (F1 + beta2)) * strain[0] + vort[0] - q[0])
                        / (m_shapeParams[1] + beta2 * m_shapeParams[3])
                  + pitch[0];
          tq[1] = q[0] * q[2] * (F1 - beta2) / (beta2 + F1)
                  + angAccFac * (((beta2 - F1) / (F1 + beta2)) * strain[1] + vort[1] - q[1])
                        / (m_shapeParams[1] + beta2 * m_shapeParams[3])
                  + pitch[1];
          tq[2] = angAccFac * (vort[2] - q[2]) / (F2 * m_shapeParams[1]) + pitch[2];
 
          for(MInt i = 0; i < 3; i++)
            rhs[i] = q[i] - q0[i] - fdtB2 * (tq0[i] + tq[i]);
 
          delta = F0;
          for(MInt i = 0; i < 3; i++) {
            const MFloat qq = q[i];
            for(MInt j = 0; j < 3; j++) {
              q[i] -= iW(i, j) * rhs(j);
            }
            delta = mMax(delta, fabs(q[i] - qq));
          }
          it++;
        }
 
        if(it >= maxit || it <= 0) {
          cerr << "Newton iterations did not converge " << m_partId << endl;
        }
 
        for(MInt i = 0; i < 3; i++) {
          m_angularVel[i] = q[i];
          m_angularAccel[i] = tq[i];
        }
 
        // integrate quaternions
        w = m_oldQuaternion[0];
        x = m_oldQuaternion[1];
        y = m_oldQuaternion[2];
        z = m_oldQuaternion[3];
 
        W(0, 0) = F0;
        W(0, 1) = -q[0];
        W(0, 2) = -q[1];
        W(0, 3) = -q[2];
        W(1, 0) = q[0];
        W(1, 1) = F0;
        W(1, 2) = q[2];
        W(1, 3) = -q[1];
        W(2, 0) = q[1];
        W(2, 1) = -q[2];
        W(2, 2) = F0;
        W(2, 3) = q[0];
        W(3, 0) = q[2];
        W(3, 1) = q[1];
        W(3, 2) = -q[0];
        W(3, 3) = F0;
 
        for(MInt i = 0; i < 4; i++) {
          for(MInt j = 0; j < 4; j++) {
            W(i, j) = -F1B2 * fdtB2 * W(i, j);
          }
          W(i, i) = F1;
        }
 
        rhs(0) = w + F1B2 * fdtB2 * (-x * q0[0] - y * q0[1] - z * q0[2]);
        rhs(1) = x + F1B2 * fdtB2 * (w * q0[0] - z * q0[1] + y * q0[2]);
        rhs(2) = y + F1B2 * fdtB2 * (z * q0[0] + w * q0[1] - x * q0[2]);
        rhs(3) = z + F1B2 * fdtB2 * (-y * q0[0] + x * q0[1] + w * q0[2]);
 
        maia::math::invert(W, iW, 4, 4);
 
        for(MInt i = 0; i < 4; i++) {
          m_quaternion[i] = F0;
          for(MInt j = 0; j < 4; j++) {
            m_quaternion[i] += iW(i, j) * rhs(j);
          }
        }
 
        MFloat abs = F0;
        for(MInt i = 0; i < 4; i++)
          abs += POW2(m_quaternion[i]);
        for(MInt i = 0; i < 4; i++)
          m_quaternion[i] /= sqrt(abs);
 
        // integrate linear motion
        K.fill(F0);
        K(0, 0) = F1 / (m_shapeParams[0] / POW2(m_semiMinorAxis) + m_shapeParams[1]);
        K(1, 1) = K(0, 0);
        K(2, 2) = F1 / (m_shapeParams[0] / POW2(m_semiMinorAxis) + beta2 * m_shapeParams[3]);
 
        // new drag correlations
        MFloatScratchSpace accDrag(3, AT_, "accDrag");
        accDrag.fill(F0);
        MFloat facCL = liftFactor(Rep, beta, inclinationAngle, K(0, 0), K(2, 2));
        MFloat facCD = dragFactor(Rep, beta, inclinationAngle, K(0, 0), K(2, 2));
        for(MInt i = 0; i < nDim; i++) {
          accLiftHat[i] = addFacCd * facCL * magRelVelHat * dirLHat[i];
          accDragHat[i] = addFacCd * facCD * magRelVelHat * dirDHat[i];
          accHat[i] = accDragHat[i] + accLiftHat[i];
        }
        maia::math::matrixVectorProductTranspose(&accDrag[0], R, &accHat[0]);
 
        for(MInt i = 0; i < nDim; i++) {
          m_accel[i] = accDrag[i] + (1.0 - invDensityRatio) * s_Frm[i];
          m_velocity[i] = m_oldVel[i] + dt * F1B2 * (m_accel[i] + m_oldAccel[i]);
          m_position[i] = m_oldPos[i] + dt * F1B2 * (m_velocity[i] + m_oldVel[i]) * s_lengthFactor
                          + F1B4 * POW2(dt) * (m_accel[i] + m_oldAccel[i]);
        }
        break;
      }
      default: {
        mTerm(1, AT_, "Unknown particle corrector equation type!");
      }
    }
  }
 
#ifdef LPT_DEBUG
  if(m_partId == debugPartId) {
    cerr << "AT " << m_velocity[0] << " " << m_velocity[1] << " " << m_velocity[2] << endl;
  }
#endif
 
  // 3. update cellId based on the corrected position and set the new status!
  checkCellChange(&m_oldPos[0]);
 
#ifdef LPT_DEBUG
  for(MInt i = 0; i < nDim; i++) {
    if(std::isnan(m_accel[i])) {
      cerr << "Nan acc. for " << s_backPtr->domainId() << " " << m_partId << " " << m_position[0] << " "
           << m_position[1] << " " << m_position[nDim - 1] << " " << s_backPtr->a_isValidCell(m_cellId) << " "
           << m_oldPos[0] << " " << m_oldPos[1] << " " << m_oldPos[nDim - 1] << " " << m_creationTime << " "
           << fluidVel[0] << " " << fluidVel[1] << " " << fluidVel[nDim - 1] << fluidDensity << " " << T << " "
           << fluidViscosity << " " << m_oldVel[0] << " " << m_oldVel[1] << " " << m_oldVel[nDim - 1] << endl;
    }
  }
#endif
}

operator=()

template<MInt nDim>

LPTEllipsoidal & LPTEllipsoidal< nDim >::operator= ( const LPTEllipsoidal< nDim > &  )

default

particleMass()

template<MInt nDim>

MFloat LPTEllipsoidal< nDim >::particleMass ( ) const

inline

Definition at line 162 of file lptellipsoidal.h.

                                     {
    return F4B3 * PI * POW3(m_semiMinorAxis) * m_aspectRatio * s_backPtr->m_material->density();
  }

particleRadius()

template<MInt nDim>

MFloat LPTEllipsoidal< nDim >::particleRadius ( MInt  direction ) const

inline

Definition at line 145 of file lptellipsoidal.h.

                                                     {
    if(direction == m_particleMajorAxis) {
      return maxParticleRadius();
    }
    return minParticleRadius();
  }

particleRe()

template<MInt nDim>

MFloat LPTEllipsoidal< nDim >::particleRe ( const MFloat  velocity, const MFloat  density, const MFloat  dynamicViscosity  )

inline

Definition at line 176 of file lptellipsoidal.h.

                                                                                                       {
    return density * velocity * m_eqDiameter / dynamicViscosity;
  }

particleVolume()

template<MInt nDim>

MFloat LPTEllipsoidal< nDim >::particleVolume ( ) const

inline

Definition at line 159 of file lptellipsoidal.h.

{ return F4B3 * PI * POW3(m_semiMinorAxis) * m_aspectRatio; }

resetWeights()

template<MInt nDim>

void LPTEllipsoidal< nDim >::resetWeights ( )

inlineoverridevirtual

Implements LPTBase< nDim >.

Definition at line 121 of file lptellipsoidal.h.

                               {
    m_neighborList.clear();
    this->template interpolateAndCalcWeights<0, 0>(m_cellId, m_position.data(), nullptr, m_redistWeight);
  };

setParticleFluidVelocities()

template<MInt nDim>

void LPTEllipsoidal< nDim >::setParticleFluidVelocities	(	MInt	cellId,
		MFloat *	position,
		MFloat *	quaternions,
		MFloat *	fluidVel,
		MFloat *	fluidVelGrad
	)

torqueFactor()

template<MInt nDim>

MFloat LPTEllipsoidal< nDim >::torqueFactor	(	const MFloat	partRe,
		const MFloat	beta,
		const MFloat	inclinationAngle
	)

Author

Laurent Andre

Date

August 2022

Definition at line 821 of file lptellipsoidal.cpp.

                                                                                                               {
  MFloat result = NAN;
  switch(s_backPtr->m_torqueModelType) {
    case 0: { // torque switched off
      result = F0;
      break;
    }
    case 1: { // linear
      result = F1;
      break;
    }
    case 2: { // new correlations by Konstantin (2020)
      std::array<MFloat, 7> tCorrs{0.931, 0.675, 0.162, 0.657, 2.77, 0.178, 0.177};
      MFloat CTMax = tCorrs[0] * pow(log(beta), tCorrs[1]) / pow(partRe, tCorrs[2])
                     + tCorrs[3] * pow(log(beta), tCorrs[4]) / pow(partRe, tCorrs[5] + tCorrs[6] * log(beta));
      MFloat CT = F2 * sin(inclinationAngle) * cos(inclinationAngle) * CTMax;
      result = CT;
      break;
    }
    default: {
      mTerm(1, AT_, "Unknown particle torque method");
    }
  }
  return result;
}

Friends And Related Function Documentation

operator<<

template<MInt nDim>

template<MInt >

std::ostream & operator<< ( std::ostream &  os, const LPTEllipsoidal< nDim > &  m  )

friend

Definition at line 243 of file lptellipsoidal.h.

                                                                      {
  return os << m.m_partId;
}

Member Data Documentation

m_angularAccel

template<MInt nDim>

std::array<MFloat, nDim> LPTEllipsoidal< nDim >::m_angularAccel {}

Definition at line 63 of file lptellipsoidal.h.

m_angularVel

template<MInt nDim>

std::array<MFloat, nDim> LPTEllipsoidal< nDim >::m_angularVel {}

Definition at line 61 of file lptellipsoidal.h.

m_aspectRatio

template<MInt nDim>

MFloat LPTEllipsoidal< nDim >::m_aspectRatio = std::numeric_limits<MFloat>::quiet_NaN()

Definition at line 85 of file lptellipsoidal.h.

m_eqDiameter

template<MInt nDim>

MFloat LPTEllipsoidal< nDim >::m_eqDiameter = std::numeric_limits<MFloat>::quiet_NaN()

private

Definition at line 209 of file lptellipsoidal.h.

m_fluidDensity

template<MInt nDim>

MFloat LPTEllipsoidal< nDim >::m_fluidDensity = std::numeric_limits<MFloat>::quiet_NaN()

Definition at line 67 of file lptellipsoidal.h.

m_fluidVel

template<MInt nDim>

std::array<MFloat, nDim> LPTEllipsoidal< nDim >::m_fluidVel {}

Definition at line 65 of file lptellipsoidal.h.

m_fluidVelMag

template<MInt nDim>

MFloat LPTEllipsoidal< nDim >::m_fluidVelMag = std::numeric_limits<MFloat>::quiet_NaN()

Definition at line 98 of file lptellipsoidal.h.

m_heatFlux

template<MInt nDim>

MFloat LPTEllipsoidal< nDim >::m_heatFlux = std::numeric_limits<MFloat>::quiet_NaN()

Definition at line 100 of file lptellipsoidal.h.

m_oldAccel

template<MInt nDim>

std::array<MFloat, nDim> LPTEllipsoidal< nDim >::m_oldAccel {}

Definition at line 70 of file lptellipsoidal.h.

m_oldAngularAccel

template<MInt nDim>

std::array<MFloat, nDim> LPTEllipsoidal< nDim >::m_oldAngularAccel {}

Definition at line 74 of file lptellipsoidal.h.

m_oldAngularVel

template<MInt nDim>

std::array<MFloat, nDim> LPTEllipsoidal< nDim >::m_oldAngularVel {}

Definition at line 72 of file lptellipsoidal.h.

m_oldFluidDensity

template<MInt nDim>

MFloat LPTEllipsoidal< nDim >::m_oldFluidDensity = std::numeric_limits<MFloat>::quiet_NaN()

Definition at line 82 of file lptellipsoidal.h.

m_oldFluidVel

template<MInt nDim>

std::array<MFloat, nDim> LPTEllipsoidal< nDim >::m_oldFluidVel {}

Definition at line 76 of file lptellipsoidal.h.

m_oldQuaternion

template<MInt nDim>

std::array<MFloat, 4> LPTEllipsoidal< nDim >::m_oldQuaternion {}

Definition at line 93 of file lptellipsoidal.h.

m_particleMajorAxis

template<MInt nDim>

MInt LPTEllipsoidal< nDim >::m_particleMajorAxis = 2

Definition at line 105 of file lptellipsoidal.h.

m_quaternion

template<MInt nDim>

std::array<MFloat, 4> LPTEllipsoidal< nDim >::m_quaternion {}

Definition at line 92 of file lptellipsoidal.h.

m_redistWeight

template<MInt nDim>

std::vector<MFloat> LPTEllipsoidal< nDim >::m_redistWeight {}

Definition at line 115 of file lptellipsoidal.h.

m_semiMinorAxis

template<MInt nDim>

MFloat LPTEllipsoidal< nDim >::m_semiMinorAxis = std::numeric_limits<MFloat>::quiet_NaN()

Definition at line 87 of file lptellipsoidal.h.

m_shapeParams

template<MInt nDim>

std::array<MFloat, 4> LPTEllipsoidal< nDim >::m_shapeParams {std::numeric_limits<MFloat>::quiet_NaN()}

Definition at line 89 of file lptellipsoidal.h.

m_temperature

template<MInt nDim>

MFloat LPTEllipsoidal< nDim >::m_temperature = std::numeric_limits<MFloat>::quiet_NaN()

Definition at line 96 of file lptellipsoidal.h.

m_velocityGradientFluid

template<MInt nDim>

std::array<MFloat, nDim * nDim> LPTEllipsoidal< nDim >::m_velocityGradientFluid {}

Definition at line 79 of file lptellipsoidal.h.

s_floatElements

template<MInt nDim>

constexpr MInt LPTEllipsoidal< nDim >::s_floatElements = 7 + 8 * nDim + 4 * 2

staticconstexpr

Definition at line 102 of file lptellipsoidal.h.

s_Frm

template<MInt nDim>

std::array< MFloat, nDim > LPTEllipsoidal< nDim >::s_Frm {}

static

Definition at line 112 of file lptellipsoidal.h.

s_intElements

template<MInt nDim>

constexpr MInt LPTEllipsoidal< nDim >::s_intElements = 0

staticconstexpr

Definition at line 103 of file lptellipsoidal.h.

s_lengthFactor

template<MInt nDim>

MFloat LPTEllipsoidal< nDim >::s_lengthFactor = F1

static

Definition at line 107 of file lptellipsoidal.h.

s_Pr

template<MInt nDim>

MFloat LPTEllipsoidal< nDim >::s_Pr = F1

static

Definition at line 109 of file lptellipsoidal.h.

s_Re

template<MInt nDim>

MFloat LPTEllipsoidal< nDim >::s_Re = F1

static

Definition at line 108 of file lptellipsoidal.h.

s_Sc

template<MInt nDim>

MFloat LPTEllipsoidal< nDim >::s_Sc = F1

static

Definition at line 110 of file lptellipsoidal.h.

s_We

template<MInt nDim>

MFloat LPTEllipsoidal< nDim >::s_We

static

Definition at line 111 of file lptellipsoidal.h.

---

The documentation for this class was generated from the following files:

• /home/gitlab-runner/scratch/builds/oxpnswJ6/1/aia/m-AIA/m-AIA/src/LPT/lpt.h
• /home/gitlab-runner/scratch/builds/oxpnswJ6/1/aia/m-AIA/m-AIA/src/LPT/lptellipsoidal.h
• /home/gitlab-runner/scratch/builds/oxpnswJ6/1/aia/m-AIA/m-AIA/src/LPT/lptellipsoidal.cpp