HeatFluidCoupledSpecies.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 HEATFLUIDCOUPLEDSPECIES_H
#define HEATFLUIDCOUPLEDSPECIES_H
#include "ThermalSpecies.h"
#include "Interactions/NormalForceInteractions/HeatFluidCoupledInteraction.h"
class BaseInteraction;
 
template<class NormalForceSpecies>
class HeatFluidCoupledSpecies : public ThermalSpecies<NormalForceSpecies>
{
public:
 
    typedef HeatFluidCoupledInteraction<typename NormalForceSpecies::InteractionType> InteractionType;
 
    HeatFluidCoupledSpecies();
 
    HeatFluidCoupledSpecies(const HeatFluidCoupledSpecies& s);
 
    virtual ~HeatFluidCoupledSpecies();
 
    void write(std::ostream& os) const;
 
    void read(std::istream& is);
 
    std::string getBaseName() const;
 
    Mdouble getMassTransferCoefficient() const;
 
    void setMassTransferCoefficient(Mdouble massTransferCoefficient);
 
    Mdouble getLatentHeatVaporization() const;
 
    void setLatentHeatVaporization(Mdouble latentHeatVaporization);
 
    Mdouble getLiquidDensity() const;
 
    void setLiquidDensity(Mdouble liquidDensity);
 
    Mdouble getEvaporationCoefficientA() const;
 
    void setEvaporationCoefficientA(Mdouble evaporationCoefficientA);
 
    Mdouble getEvaporationCoefficientB() const;
 
    void setEvaporationCoefficientB(Mdouble evaporationCoefficientB);
 
    Mdouble getAmbientHumidity() const;
 
    void setAmbientHumidity(Mdouble ambientHumidity);
 
    Mdouble getAmbientEquilibriumMoistureContent() const;
 
    void setAmbientEquilibriumMoistureContent(Mdouble ambientEquilibriumMoistureContent);
 
    Mdouble getAmbientVapourConcentration() const;
 
    void setAmbientVapourConcentration(Mdouble ambientVapourConcentration);
 
    void actionsAfterTimeStep(BaseParticle* particle) const override;
 
    std::array<double,2> f(double liquidVolume,double temperature, double mass, double surfaceArea) const;
 
private:
 
    Mdouble massTransferCoefficient_;
 
    Mdouble latentHeatVaporization_;
 
    Mdouble liquidDensity_;
 
    Mdouble evaporationCoefficientA_;
 
    Mdouble evaporationCoefficientB_;
 
    Mdouble ambientHumidity_;
 
    Mdouble ambientEquilibriumMoistureContent_;
 
    Mdouble ambientVapourConcentration_;
 
};
 
template<class NormalForceSpecies>
HeatFluidCoupledSpecies<NormalForceSpecies>::HeatFluidCoupledSpecies()
        : ThermalSpecies<NormalForceSpecies>()
{
    massTransferCoefficient_=0.0;
    latentHeatVaporization_=0.0;
    liquidDensity_=1.0;
    evaporationCoefficientA_=0.0;
    evaporationCoefficientB_=0.0;
    // note, this default value is unrealistically high; 0.3 is realistic.
    ambientHumidity_=1.0;
    ambientEquilibriumMoistureContent_=0.0;
    ambientVapourConcentration_=0.0;
}
 
template<class NormalForceSpecies>
HeatFluidCoupledSpecies<NormalForceSpecies>::HeatFluidCoupledSpecies(const HeatFluidCoupledSpecies& s)
        : ThermalSpecies<NormalForceSpecies>(s)
{
    massTransferCoefficient_= s.massTransferCoefficient_;
    latentHeatVaporization_= s.latentHeatVaporization_;
    liquidDensity_=s.liquidDensity_;
    evaporationCoefficientA_=s.evaporationCoefficientA_;
    evaporationCoefficientB_=s.evaporationCoefficientB_;
    ambientHumidity_=s.ambientHumidity_;
    ambientEquilibriumMoistureContent_=s.ambientEquilibriumMoistureContent_;
    ambientVapourConcentration_=s.ambientVapourConcentration_;
}
 
template<class NormalForceSpecies>
HeatFluidCoupledSpecies<NormalForceSpecies>::~HeatFluidCoupledSpecies()
{}
 
template<class NormalForceSpecies>
void HeatFluidCoupledSpecies<NormalForceSpecies>::write(std::ostream& os) const
{
    ThermalSpecies<NormalForceSpecies>::write(os);
    os << " massTransferCoefficient " << massTransferCoefficient_;
    os << " latentHeatVaporization " << latentHeatVaporization_;
    os << " liquidDensity " << liquidDensity_;
    os << " evaporationCoefficientA " << evaporationCoefficientA_;
    os << " evaporationCoefficientA " << evaporationCoefficientB_;
    os << " ambientHumidity " << ambientHumidity_;
    os << " ambientEquilibriumMoistureContent " << ambientEquilibriumMoistureContent_;
    os << " ambientVapourConcentration " << ambientVapourConcentration_;
}
 
template<class NormalForceSpecies>
void HeatFluidCoupledSpecies<NormalForceSpecies>::read(std::istream& is)
{
    std::string dummy;
    ThermalSpecies<NormalForceSpecies>::read(is);
    is >> dummy >> massTransferCoefficient_;
    is >> dummy >> latentHeatVaporization_;
    is >> dummy >> liquidDensity_;
    is >> dummy >> evaporationCoefficientA_;
    is >> dummy >> evaporationCoefficientB_;
    is >> dummy >> ambientHumidity_;
    is >> dummy >> ambientEquilibriumMoistureContent_;
    is >> dummy >> ambientVapourConcentration_;
    double ambientTemperature;
    if (helpers::readOptionalVariable<Mdouble>(is, "ambientTemperature", ambientTemperature)) {
        ThermalSpecies<NormalForceSpecies>::setAmbientTemperature(ambientTemperature);
    }
}
 
template<class NormalForceSpecies>
std::string HeatFluidCoupledSpecies<NormalForceSpecies>::getBaseName() const
{
    return "HeatFluidCoupled" + NormalForceSpecies::getBaseName();
}
 
template<class NormalForceSpecies>
Mdouble HeatFluidCoupledSpecies<NormalForceSpecies>::getMassTransferCoefficient() const
{
    return massTransferCoefficient_;
}
 
template<class NormalForceSpecies>
void HeatFluidCoupledSpecies<NormalForceSpecies>::setMassTransferCoefficient(Mdouble massTransferCoefficient)
{
    logger.assert_always(massTransferCoefficient > 0,
                         "[HeatFluidCoupledSpecies<>::setMassTransferCoefficient(%)] value has to be positive",
                         massTransferCoefficient);
    massTransferCoefficient_ = massTransferCoefficient;
}
 
template<class NormalForceSpecies>
Mdouble HeatFluidCoupledSpecies<NormalForceSpecies>::getLatentHeatVaporization() const
{
    return latentHeatVaporization_;
}
 
template<class NormalForceSpecies>
void HeatFluidCoupledSpecies<NormalForceSpecies>::setLatentHeatVaporization(Mdouble latentHeatVaporization)
{
    logger.assert_always(latentHeatVaporization > 0,
                         "[HeatFluidCoupledSpecies<>::setLatentHeatVaporization(%)] value has to be positive",
                         latentHeatVaporization);
    latentHeatVaporization_ = latentHeatVaporization;
}
 
template<class NormalForceSpecies>
Mdouble HeatFluidCoupledSpecies<NormalForceSpecies>::getLiquidDensity() const
{
    return liquidDensity_;
}
 
template<class NormalForceSpecies>
void HeatFluidCoupledSpecies<NormalForceSpecies>::setLiquidDensity(Mdouble liquidDensity)
{
    logger.assert_always(liquidDensity > 0,
                         "[HeatFluidCoupledSpecies<>::setLiquidDensity(%)] value has to be positive",
                         liquidDensity);
    liquidDensity_ = liquidDensity;
}
 
template<class NormalForceSpecies>
Mdouble HeatFluidCoupledSpecies<NormalForceSpecies>::getEvaporationCoefficientA() const
{
    return evaporationCoefficientA_;
}
 
template<class NormalForceSpecies>
void HeatFluidCoupledSpecies<NormalForceSpecies>::setEvaporationCoefficientA(Mdouble evaporationCoefficientA)
{
    logger.assert_always(evaporationCoefficientA >= 0,
                         "[HeatFluidCoupledSpecies<>::setEvaporationCoefficientA(%)] value has to be positive",
                         evaporationCoefficientA);
    evaporationCoefficientA_ = evaporationCoefficientA;
}
 
template<class NormalForceSpecies>
Mdouble HeatFluidCoupledSpecies<NormalForceSpecies>::getEvaporationCoefficientB() const
{
    return evaporationCoefficientB_;
}
 
template<class NormalForceSpecies>
void HeatFluidCoupledSpecies<NormalForceSpecies>::setEvaporationCoefficientB(Mdouble evaporationCoefficientB)
{
    logger.assert_always(evaporationCoefficientB <= 0,
                         "[HeatFluidCoupledSpecies<>::setEvaporationCoefficientB(%)] value has to be negative",
                         evaporationCoefficientB);
    evaporationCoefficientB_ = evaporationCoefficientB;
}
 
template<class NormalForceSpecies>
Mdouble HeatFluidCoupledSpecies<NormalForceSpecies>::getAmbientHumidity() const
{
    return ambientHumidity_;
}
 
template<class NormalForceSpecies>
void HeatFluidCoupledSpecies<NormalForceSpecies>::setAmbientHumidity(Mdouble ambientHumidity)
{
    //note, this is not allowed to be zero!
    logger.assert_always(ambientHumidity > 0,
                         "[HeatFluidCoupledSpecies<>::setAmbientHumidity(%)] value has to be positive",
                         ambientHumidity);
    ambientHumidity_ = ambientHumidity;
}
 
template<class NormalForceSpecies>
Mdouble HeatFluidCoupledSpecies<NormalForceSpecies>::getAmbientEquilibriumMoistureContent() const
{
    return ambientEquilibriumMoistureContent_;
}
 
template<class NormalForceSpecies>
void HeatFluidCoupledSpecies<NormalForceSpecies>::setAmbientEquilibriumMoistureContent(Mdouble ambientEquilibriumMoistureContent)
{
    logger.assert_always(ambientEquilibriumMoistureContent >= 0,
                         "[HeatFluidCoupledSpecies<>::setAmbientEquilibriumMoistureContent(%)] value has to be positive",
                         ambientEquilibriumMoistureContent);
    ambientEquilibriumMoistureContent_ = ambientEquilibriumMoistureContent;
}
 
template<class NormalForceSpecies>
Mdouble HeatFluidCoupledSpecies<NormalForceSpecies>::getAmbientVapourConcentration() const
{
    return ambientVapourConcentration_;
}
 
template<class NormalForceSpecies>
void HeatFluidCoupledSpecies<NormalForceSpecies>::setAmbientVapourConcentration(Mdouble ambientVapourConcentration)
{
    logger.assert_always(ambientVapourConcentration >= 0,
                         "[HeatFluidCoupledSpecies<>::setAmbientVapourConcentration(%)] value has to be positive",
                         ambientVapourConcentration);
    ambientVapourConcentration_ = ambientVapourConcentration;
}
 
template<class NormalForceSpecies>
void HeatFluidCoupledSpecies<NormalForceSpecies>::actionsAfterTimeStep(BaseParticle* baseParticle) const
{
    double dt = baseParticle->getHandler()->getDPMBase()->getTimeStep();
    auto p = dynamic_cast<HeatFluidCoupledParticle*>(baseParticle);
    double mass = p->getMass();
    double surfaceArea = p->getSurfaceArea();
 
    // Runge–Kutta method.
    std::array<double,2> k1 = f(p->getLiquidVolume()               , p->getTemperature()               , mass, surfaceArea);
    std::array<double,2> k2 = f(p->getLiquidVolume() + 0.5*dt*k1[0], p->getTemperature() + 0.5*dt*k1[1], mass, surfaceArea);
    std::array<double,2> k3 = f(p->getLiquidVolume() + 0.5*dt*k2[0], p->getTemperature() + 0.5*dt*k2[1], mass, surfaceArea);
    std::array<double,2> k4 = f(p->getLiquidVolume() +     dt*k3[0], p->getTemperature() +     dt*k3[1], mass, surfaceArea);
    double dLiquidVolume = dt * (k1[0] + 2*k2[0] + 2*k3[0] + k4[0]) / 6.0;
    double dTemperature  = dt * (k1[1] + 2*k2[1] + 2*k3[1] + k4[1]) / 6.0;
    // Change sign, so we can think in terms of subtracting a positive number, instead of adding a negative number.
    // Keep in mind, it is still possible that we are adding liquid/temperature (so now subtracting a negative number),
    // however, since there is only a lower limit of 0 but no upper limit, this is easier to think about.
    dLiquidVolume = -dLiquidVolume;
    dTemperature = -dTemperature;
 
    // Subtract as much as possible from the liquid film. Any remainder is subtracted from the liquid bridges.
    double dLiquidFilm = std::min(dLiquidVolume, p->getLiquidVolume());
    double remainder = dLiquidVolume - dLiquidFilm;
 
    if (remainder > 0.0)
    {
        // Get the total liquid bridge volume. Note, don't use p->getLiquidBridgeVolume() as it halves the bridges with
        // other particles. The order of handling the particles matters for how much liquid bridge is still available
        // to them. It might seem more correct to use half the bridge instead of the full bridge, however that would
        // only add an error, because the second particle halving the bridge would include the evaporated volume of the
        // first particle. So instead of using L/2 it would use (L-evap)/2, which completely undoes the "correction" of
        // halving the bridges in the first place.
        double liquidBridgeVolume = 0.0;
        for (auto i : p->getInteractions())
        {
            auto j = dynamic_cast<LiquidMigrationWilletInteraction*>(i);
            liquidBridgeVolume += j->getLiquidBridgeVolume();
        }
        double dLiquidBridge = std::min(remainder, liquidBridgeVolume);
        remainder -= dLiquidBridge;
 
        double factor = liquidBridgeVolume == 0.0 ? 0.0 : dLiquidBridge / liquidBridgeVolume;
        for (auto i : p->getInteractions())
        {
            auto j = dynamic_cast<LiquidMigrationWilletInteraction*>(i);
            j->addLiquidBridgeVolumeByEvaporation(-factor * j->getLiquidBridgeVolume());
        }
    }
 
    p->addLiquidVolume(-dLiquidFilm);
    p->addTotalEvaporatedLiquidVolume(dLiquidVolume - remainder);
    // If not all liquid could be removed, correct the temperature to be removed.
    double dLiquidVolumeFactor = dLiquidVolume == 0.0 ? 0.0 : (dLiquidVolume - remainder) / dLiquidVolume;
    p->setTemperature(std::max(0.0, p->getTemperature() - dTemperature * dLiquidVolumeFactor));
}
 
template<class NormalForceSpecies>
std::array<double, 2> HeatFluidCoupledSpecies<NormalForceSpecies>::f(double liquidVolume, double temperature, double mass, double surfaceArea) const
{
    // Saturated vapour concentration  (kg/m^3), eq(7)
    // Dependency on temperature is valid between 1 and 200 degree Celsius
    // (it is unphysically decreasing between 0 and 1 degree),
    // extrapolated as linear functions
    double dT = temperature - 273;
    double saturatedVapourConcentration = dT < 1 ? 8.319815774e-3 * temperature / 274 :
            (((4.844e-9*dT - 1.4807e-7)*dT + 2.6572e-5)*dT - 4.8613e-5)*dT + 8.342e-3;
 
    // Relative activation energy of evaporation, eq(6)
    double equilibriumActivationEnergy = -constants::R * this->getAmbientTemperature() * log(getAmbientHumidity());
    // \todo what is the difference between the variable X and Y in Azmir?
    double moistureContent = liquidVolume * getLiquidDensity() / mass;
    double activationEnergy = equilibriumActivationEnergy * (moistureContent - getAmbientEquilibriumMoistureContent());
 
    // Vapour concentrations at the particle-medium interface (kg/m^3), eq(5)
    // Implemented such that the code does not necessarily fail if temperature is 0
    double interfaceVapourConcentration = temperature == 0 ? 0.0 :
            exp(-activationEnergy / (constants::R * temperature)) * saturatedVapourConcentration;
 
    // Drying rate of liquid film mass (kg/s), eq (4)
    double dLiquidMass = -getMassTransferCoefficient() * surfaceArea * (interfaceVapourConcentration - getAmbientVapourConcentration());
 
    // Drying rate of liquid film volume (kg/s), eq (3)
    double dLiquidVolume = dLiquidMass / getLiquidDensity();
 
    // Heat of evaporation (J/s), eq (2)
    double heatOfEvaporation = getLatentHeatVaporization() * (1.0 + getEvaporationCoefficientA() *
            exp((getEvaporationCoefficientB() * getLiquidDensity() / mass) * liquidVolume)) * dLiquidMass;
 
    // Cooling rate of particle (K/s), eq (1)
    double dTemperature = heatOfEvaporation / (mass * this->getHeatCapacity());
 
    return { dLiquidVolume, dTemperature };
}
#endif