hydram.h

Go to the documentation of this file.

/* TODO:
 *    4/2/2004
 *     need to return a matrix of derivatives
 *      dHeat/dT, dHeat/dw
 *      dWater/dT, dWater/dw
 *     (or something like that) to enable quadratic convergence of the global iteration with respect to the heat and water internal source.
 *    now I have only double dHeatdT(ksi, T, dt) - internal, no water dependence. To master material, the hydration model must provide the derivatives based on the internal state of hydration.
 *
 *    12/1:
 *    - track versions - save intermediate versions of modified files in separate folders!
 *    11/1:
 *    - where to choose, which HydrationModel to create? HydrationModelInterface->constructor can't know the type of hydram
 *      but HydrationModelInterface->initializeFrom(ir) can select and create the appropriate hydration model
 *      maybe it should be added to usrdefsub.C
 *    x try to eliminate HydrationModelStatusInterface, to enable use of different status classes without redefining (sm / tm)
 *      will not do, master status must have reference to HydrationModelStatus
 *      tm status will be defined by simply inheriting from TransportMaterialStatus && HydrationModelStatusInterface
 *    - where to put gp hydration degree access routines? Master material->hmInterface->giveTempKsi or hydrationModel->giveTempHydrationDegree(gp)
 *    10/1:
 *    - cleanly: base class HydrationModelStatus + HydrationModel, defining interface data and services, derive HellmichHydrationModel etc.
 *    - water content influence on hydration -
 *      moisture consumption (stoichiometry of c-s-h?) -> computeInternalSourceVector(+heat,-water content)
 *    - affinity function from CEMHYD3D results: ksi(t)-> dksidt(ksi)->A(ksi) - \
 *    - real-time from CEMHYD3D (better not all IP, only selected range)      -/ needs CEMHYD3D input preprocessor for generating microstructure (VCCTL?) or input microstructure data
 *    - or add alcali reaction
 *    - hydration interpolated from CEMHYD3D with moisture and temperature
 *      needs several hydration model / master material managed statuses and services for interpolating hydration degree for given coordinates
 */
 
/*
 * USAGE:
 *
 * use multiple inheritance from HydrationModelInterface and HydrationModelStatusInterface.
 *
 * add HydrationModelInterface to material:
 * constructor
 * initializeFrom
 * saveContext / restoreContext
 * TransportMaterial::giveCharacteristicValue
 * use HMI::updateInternalState(tempStateVector, GaussPoint, TimeStep) to update the hydration status
 *
 * add HydrationModelStatusInterface to status:
 * constructor
 * updateYourself
 * printOutputAt
 * add giveInterface
 *
 * then the following services are available:
 * giveHydrationDegree(GaussPoint, TimeStep, ValueModeType)
 * computeInternalSourceVector(sourceVector, GaussPoint, TimeStep, ValueModeType)
 * computeIntSourceLHSCoeff(MAtResponseMode, GaussPoint, TimeStep)
 *
 * missing - consistent tangent, see TODO
 *
 */
 
// class HydrationModel - returns hydration degree ... dksi(ksi, dt, T) ... localResidual, Affinity
// ( - for non-isothermal analysis               hydrationHeat, capacitycoeff(ksi, dt, T) ... dksidT, dksidh, dAdksi
 
// Concrete hydration model - Hellmich
// Material extension, derived from Material to use MateriaStatus services
 
// material->computeInternalSourceVector(val, gp, tStep, mode)
//                                       answer,     chartypemode - total, velocity, acceleration, incremental
// = le*dksi
 
// capacity ... material->giveCharacteristicValue(rmode, gp, tStep)
//              rmode: material responsemode = capacity
// = rc - le*dksidT
 
 
#ifndef hydram_h
#define hydram_h
 
#include "floatarray.h"
#include "floatmatrix.h"
#include "timestep.h"
#include "material.h"
#include "interface.h"
#include "valuemodetype.h"
 
#include <memory>


#define _IFT_HydrationModel_Name "hydrationmodel" 
#define _IFT_HydrationModel_hydration "hydration"
#define _IFT_HydrationModel_c60mix "c60mix"
#define _IFT_HydrationModel_timeScale "timescale"
#define _IFT_HydrationModel_hheat "hheat"
#define _IFT_HydrationModel_cv "cv"
#define _IFT_HydrationModel_water "water"


#define _IFT_HydrationModelInterface_hydration "hydration"
#define _IFT_HydrationModelInterface_castAt "castat"
 
namespace oofem {
#define ROOT_PRECISION_DKSI 1e-14
#define BINARY_TREE_STEPS 2
 
// default time steps for evaluating the hydration degree increment
#define HYDRATION_MAXSTEP0 3600 //600
#define HYDRATION_MAXSTEP1 86400 //3600
 
 
// =========== Hydration Model Status class ============
class HydrationModelStatus : public MaterialStatus
{
protected:
    // hydration degree at beginning of current time step
    double hydrationDegree = 0.;
    // hydration degree at end of current time step (or during equilibrium iteration)
    double tempHydrationDegree = 0.;
 
public:
    HydrationModelStatus(GaussPoint * g);

    double giveTempHydrationDegree() const { return tempHydrationDegree; }
    double giveHydrationDegree() const { return hydrationDegree; }
    void setHydrationDegree(double v) { hydrationDegree = v; }
    void setTempHydrationDegree(double v) { tempHydrationDegree = v; }
 
    void initTempStatus() override;
    void updateYourself(TimeStep *tStep) override;
 
    void printOutputAt(FILE *file, TimeStep *tStep) const override;
    void saveContext(DataStream &stream, ContextMode mode) override;
    void restoreContext(DataStream &stream, ContextMode mode) override;
 
    const char *giveClassName() const override { return "HydrationModelStatus"; }
};
 
 
// =========== Hydration Model class ============
 
enum FindRootMethod { frRegula = 1, frBinTree = 2, frMixed = 3 };
enum MixtureType { mtLafarge = 1, mtHuber = 2, mtC60, mtC100 };
 
// derived from material to use status services
class HydrationModel : public Material
{
protected:
    MixtureType mixture;
    double hydrationStartMaxStep = 0., hydrationEndMaxStep = 0.;
    double initialHydrationDegree = 0.;
    double timeScale = 0.;
 
    // === Material parameters ===
    double aa = 0., 
           ba = 0.,
           ca = 0.,
           da = 0.,
 
           e0 = 0., 
           ear = 0., 
           le = 0., 
           cv = 0., 
           we = 0.; 
 
    // === Hydration degree increment evaluation ===
    // auxiliary values to enable external root finding method without passing them as parameters in each call
    double auxksi, auxdt, auxT, auxh;
    double localResidual(double dks) const; // G(dksi) 4.12
 
    // Auxiliary functions - root finding
    double regulafindroot() const;
    double bintreefindroot() const;
    double mixedfindroot() const;
 
    // === Material functions ===
    double affinity(double ksi) const;
    double dAdksi(double ksi) const;
    double dksidT(double ksi, double T, double h, double dt) const;
    double dksidh(double ksi, double T, double h, double dt) const;

    double computeIntSource(double T, double h, GaussPoint *gp, TimeStep *tStep, MatResponseMode rmode) const;
    double computeHydrationDegreeIncrement(double ksi, double T, double h, double dt);
 
public:
    FindRootMethod useFindRoot;
    HydrationModel();
    HydrationModel(MixtureType mix, FindRootMethod usefr);
 
    // === Hydration model services ===

    void setMixture(MixtureType mix);

    void initializeFrom(InputRecord &ir) override;

    double giveHydrationDegree(GaussPoint *gp, TimeStep *tStep, ValueModeType mode) const;
 
    // how to determine whether hydration degree is uptodate?
    // - in tm, this can be ensured by calling computeHydrationDegreeIncrement when updateInternalState is used
    // - in sm, the tm state changes only at the start of each step (now updateAuxState - 'auxiliary', because the data is copied from tm analysis)
    // where to get temperature / moisture? from gp->status->giveTempStateVector(val) => val[T,w]
    //   => implement updateInternalState here, that is called from master material updateInternalState/updateAuxState
    // use giveTempStateVector in HellmichMaterial
    virtual void updateInternalState(const FloatArray &vec, GaussPoint *gp, TimeStep *tStep);

    MaterialStatus *giveStatus(GaussPoint *gp) const override;

    void computeInternalSourceVector(FloatArray &val, GaussPoint *gp, TimeStep *tStep, ValueModeType mode) const;
    // } end new 5.1.2004
    virtual double _giveCharacteristicValue(double T, double h, MatResponseMode rmode, GaussPoint *gp, TimeStep *tStep) const;
    // --- identification and auxiliary functions ---
    const char *giveInputRecordName() const override { return _IFT_HydrationModel_Name; }
    const char *giveClassName() const override { return "HydrationModel"; }

    std::unique_ptr<MaterialStatus> CreateStatus(GaussPoint *gp) const override;
};
 
// =========== Interfaces for materials using the hydration model ============
class HydrationModelStatusInterface : public Interface
{
protected:
    std :: unique_ptr< MaterialStatus > hydrationModelStatus;
 
public:
    HydrationModelStatusInterface() {}

    HydrationModelStatus *giveHydrationModelStatus() { return static_cast<HydrationModelStatus*> (hydrationModelStatus.get()); }
    HydrationModelStatus* setHydrationModelStatus(std::unique_ptr<MaterialStatus> s) { hydrationModelStatus = std::move(s); return static_cast<HydrationModelStatus*> (hydrationModelStatus.get()); }

    void updateYourself(TimeStep *tStep);
    void printOutputAt(FILE *file, TimeStep *tStep) const;
};
 
class HydrationModelInterface : public Interface
{
    // interface for material access to HydrationModel
protected:
    std :: unique_ptr< HydrationModel > hydrationModel;
    double castAt = 0.;
    double constantHydrationDegree = 0.;
 
public:
    HydrationModel *giveHydrationModel() { return hydrationModel.get(); }
    virtual void initializeFrom(InputRecord &ir);
 
    void saveContext(DataStream &stream, ContextMode mode)
    {
        if ( hydrationModel ) {
            hydrationModel->saveContext(stream, mode);
        }
    }
    void restoreContext(DataStream &stream, ContextMode mode)
    {
        if ( hydrationModel ) {
            hydrationModel->restoreContext(stream, mode);
        }
    }

    virtual void updateInternalState(const FloatArray &vec, GaussPoint *gp, TimeStep *tStep);
    double giveHydrationDegree(GaussPoint *gp, TimeStep *tStep, ValueModeType mode) const;
};
} // end namespace oofem
#endif // hydram_h