ThermalSpecies.h

Go to the documentation of this file.

// This file is part of the MercuryDPM project (https://www.mercurydpm.org).
// Copyright (c), The MercuryDPM Developers Team. All rights reserved.
// License: BSD 3-Clause License; see the LICENSE file in the root directory.
 
#ifndef THERMALSPECIES_H
#define THERMALSPECIES_H
 
#include "Interactions/NormalForceInteractions/ThermalInteraction.h"
 
class BaseInteraction;
 
template<class NormalForceSpecies>
class ThermalSpecies : public NormalForceSpecies
{
public:
    
    typedef ThermalInteraction<typename NormalForceSpecies::InteractionType> InteractionType;
    
    ThermalSpecies();
    
    ThermalSpecies(const ThermalSpecies& s);
    
    virtual ~ThermalSpecies();
    
    void write(std::ostream& os) const;
    
    void read(std::istream& is);
    
    std::string getBaseName() const;
 
    void actionsAfterTimeStep(BaseParticle* particle) const override;
 
    Mdouble getHeatCapacity() const;
    
    void setHeatCapacity(Mdouble heatCapacity);
    
    Mdouble getThermalConductivity() const;
    
    void setThermalConductivity(Mdouble thermalConductivity);
 
    Mdouble getAmbientTemperature() const
    {
        return ambientTemperature_;
    }
 
    void setAmbientTemperature(Mdouble ambientTemperature)
    {
        logger.assert_always(ambientTemperature >= 0,"ambientTemperature cannot be negative");
        ambientTemperature_ = ambientTemperature;
    }
 
    Vec3D getAmbientVelocity() const {
        return ambientVelocity_;
    }
 
    void setAmbientVelocity(Vec3D ambientVelocity) {
        ambientVelocity_ = ambientVelocity;
    }
 
    Mdouble getAmbientThermalConductivity() const {
        return ambientThermalConductivity_;
    }
 
    void setAmbientThermalConductivity(Mdouble ambientThermalConductivity) {
        logger.assert_always(ambientThermalConductivity>0,"ambientThermalConductivity has to be positive");
        ambientThermalConductivity_ = ambientThermalConductivity;
    }
 
    Mdouble getAmbientDynamicViscosity() const {
        return ambientDynamicViscosity_;
    }
 
    void setAmbientDynamicViscosity(Mdouble ambientDynamicViscosity) {
        logger.assert_always(ambientDynamicViscosity>0,"ambientDynamicViscosity has to be positive");
        ambientDynamicViscosity_ = ambientDynamicViscosity;
    }
 
    Mdouble getAmbientHeatCapacity() const {
        return ambientHeatCapacity_;
    }
 
    void setAmbientHeatCapacity(Mdouble ambientHeatCapacity) {
        logger.assert_always(ambientHeatCapacity>0,"ambientHeatCapacity has to be positive");
        ambientHeatCapacity_ = ambientHeatCapacity;
    }
 
    Mdouble getAmbientDensity() const {
        return ambientDensity_;
    }
 
    void setAmbientDensity(Mdouble ambientDensity) {
        logger.assert_always(ambientDensity>0,"ambientDensity has to be positive");
        ambientDensity_ = ambientDensity;
    }
 
 
private:
    Mdouble heatCapacity_;
    
    Mdouble thermalConductivity_;
 
    Mdouble ambientTemperature_;
 
    // heat transfer properties
 
    Vec3D ambientVelocity_;
 
    double ambientThermalConductivity_;
 
    double ambientDynamicViscosity_;
 
    double ambientHeatCapacity_;
 
    double ambientDensity_;
};
 
template<class NormalForceSpecies>
ThermalSpecies<NormalForceSpecies>::ThermalSpecies()
        : NormalForceSpecies()
{
    heatCapacity_ = 0.0;
    thermalConductivity_ = 0.0;
    ambientTemperature_=0.0;
    ambientVelocity_ = {0,0,0};
    ambientThermalConductivity_ = 0.0;
    ambientDynamicViscosity_ = 0.0;
    ambientHeatCapacity_ = 0.0;
    ambientDensity_ = 0.0;
}
 
template<class NormalForceSpecies>
ThermalSpecies<NormalForceSpecies>::ThermalSpecies(const ThermalSpecies& s)
        : NormalForceSpecies(s)
{
    heatCapacity_ = s.heatCapacity_;
    thermalConductivity_ = s.thermalConductivity_;
    ambientTemperature_=s.ambientTemperature_;
    ambientVelocity_ = s.ambientVelocity_;
    ambientThermalConductivity_ = s.ambientThermalConductivity_;
    ambientDynamicViscosity_ = s.ambientDynamicViscosity_;
    ambientHeatCapacity_ = s.ambientHeatCapacity_;
    ambientDensity_ = s.ambientDensity_;
}
 
template<class NormalForceSpecies>
ThermalSpecies<NormalForceSpecies>::~ThermalSpecies()
{}
 
template<class NormalForceSpecies>
void ThermalSpecies<NormalForceSpecies>::write(std::ostream& os) const
{
    NormalForceSpecies::write(os);
    os << " heatCapacity " << heatCapacity_;
    os << " thermalConductivity " << thermalConductivity_;
    os << " ambientTemperature " << ambientTemperature_;
    os << " ambientVelocity " << ambientVelocity_;
    os << " ambientThermalConductivity " << ambientThermalConductivity_;
    os << " ambientDynamicViscosity " << ambientDynamicViscosity_;
    os << " ambientHeatCapacity " << ambientHeatCapacity_;
    os << " ambientDensity " << ambientDensity_;
}
 
template<class NormalForceSpecies>
void ThermalSpecies<NormalForceSpecies>::read(std::istream& is)
{
    std::string dummy;
    NormalForceSpecies::read(is);
    is >> dummy >> heatCapacity_;
    is >> dummy >> thermalConductivity_;
    helpers::readOptionalVariable<Mdouble>(is, "ambientTemperature", ambientTemperature_);
    helpers::readOptionalVariable<Vec3D>(is, "ambientVelocity", ambientVelocity_);
    helpers::readOptionalVariable<Mdouble>(is, "ambientThermalConductivity", ambientThermalConductivity_);
    helpers::readOptionalVariable<Mdouble>(is, "ambientDynamicViscosity", ambientDynamicViscosity_);
    helpers::readOptionalVariable<Mdouble>(is, "ambientHeatCapacity", ambientHeatCapacity_);
    helpers::readOptionalVariable<Mdouble>(is, "ambientDensity", ambientDensity_);
}
 
template<class NormalForceSpecies>
std::string ThermalSpecies<NormalForceSpecies>::getBaseName() const
{
    return "Thermal" + NormalForceSpecies::getBaseName();
}
 
 
template<class NormalForceSpecies>
Mdouble ThermalSpecies<NormalForceSpecies>::getHeatCapacity() const
{
    return heatCapacity_;
}
 
template<class NormalForceSpecies>
void ThermalSpecies<NormalForceSpecies>::setHeatCapacity(Mdouble heatCapacity)
{
    logger.assert_always(heatCapacity > 0,
                         "[ThermalSpecies<>::setHeatCapacity(%)] value has to be positive",
                         heatCapacity);
    heatCapacity_ = heatCapacity;
}
 
template<class NormalForceSpecies>
Mdouble ThermalSpecies<NormalForceSpecies>::getThermalConductivity() const
{
    return thermalConductivity_;
}
 
template<class NormalForceSpecies>
void ThermalSpecies<NormalForceSpecies>::setThermalConductivity(Mdouble thermalConductivity)
{
    logger.assert_always(thermalConductivity >= 0,
                         "[ThermalSpecies<>::setThermalConductivity(%)] value has to be positive",
                         thermalConductivity);
    thermalConductivity_ = thermalConductivity;
}
 
template<class NormalForceSpecies>
void ThermalSpecies<NormalForceSpecies>::actionsAfterTimeStep(BaseParticle* baseParticle) const {
    /* Compute temperature change due to heat convection with the ambient fluid */
    if (getAmbientThermalConductivity()>0)
    {
        auto p = dynamic_cast<ThermalParticle*>(baseParticle);
        double diameter = 2.0 * p->getRadius();
        double prandtl = getAmbientDynamicViscosity() * getAmbientHeatCapacity() /
                         getAmbientThermalConductivity();
        double reynolds = getAmbientDensity() * (p->getVelocity() - getAmbientVelocity()).getLength() * diameter / getAmbientDynamicViscosity();
        double nusselt = 2.0 + 0.4 * sqrt(reynolds) * cbrt(prandtl);
        double heatTransferCoefficient = nusselt * getAmbientThermalConductivity() / diameter;
        double heatTransferRate = heatTransferCoefficient * p->getSurfaceArea() *
                                  (getAmbientTemperature()-p->getTemperature());
        Mdouble temperatureChange = heatTransferRate * p->getInvMass() / getHeatCapacity()
                                    * p->getHandler()->getDPMBase()->getTimeStep();
        p->addTemperature(temperatureChange);
    }
}
 
#endif