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;
    // hydration degree at end of current time step (or during equilibrium iteration)
    double tempHydrationDegree;

public:
    HydrationModelStatus(int n, Domain * d, GaussPoint * g);
    virtual ~HydrationModelStatus() { }

    double giveTempHydrationDegree() { return tempHydrationDegree; }
    double giveHydrationDegree() { return hydrationDegree; }
    void setHydrationDegree(double v) { hydrationDegree = v; }
    void setTempHydrationDegree(double v) { tempHydrationDegree = v; }

    virtual void initTempStatus();
    virtual void updateYourself(TimeStep *tStep);

    virtual void printOutputAt(FILE *file, TimeStep *tStep);
    virtual contextIOResultType saveContext(DataStream &stream, ContextMode mode, void *obj = NULL);
    virtual contextIOResultType restoreContext(DataStream &stream, ContextMode mode, void *obj = NULL);

    // --- identification and auxiliary functions ---
    virtual const char *giveClassName() const { 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, hydrationEndMaxStep;
    double initialHydrationDegree;
    double timeScale;

    // === Material parameters ===
    double aa, 
           ba,
           ca,
           da,

           e0, 
           ear, 
           le, 
           cv, 
           we; 

    // === 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); // G(dksi) 4.12

    // Auxiliary functions - root finding
    double regulafindroot();
    double bintreefindroot();
    double mixedfindroot();

    // === Material functions ===
    double affinity(double ksi);
    double dAdksi(double ksi);
    double dksidT(double ksi, double T, double h, double dt);
    double dksidh(double ksi, double T, double h, double dt);

    double computeIntSource(const FloatArray &vec, GaussPoint *gp, TimeStep *tStep, MatResponseMode rmode);
    double computeHydrationDegreeIncrement(double ksi, double T, double h, double dt);

public:
    FindRootMethod useFindRoot;
    HydrationModel();
    HydrationModel(MixtureType mix, FindRootMethod usefr);
    virtual ~HydrationModel() { }


    // === Hydration model services ===

    void setMixture(MixtureType mix);

    virtual IRResultType initializeFrom(InputRecord *ir);

    double giveHydrationDegree(GaussPoint *gp, TimeStep *tStep, ValueModeType mode);

    // 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);

    virtual MaterialStatus *giveStatus(GaussPoint *gp) const;

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

protected:
    virtual MaterialStatus *CreateStatus(GaussPoint *gp) const;
};

// =========== Interfaces for materials using the hydration model ============
class HydrationModelStatusInterface : public Interface
{
protected:
    std :: unique_ptr< HydrationModelStatus > hydrationModelStatus;
public:
    HydrationModelStatusInterface() {}
    virtual ~HydrationModelStatusInterface() {}

    HydrationModelStatus *giveHydrationModelStatus() { return hydrationModelStatus.get(); }
    void setHydrationModelStatus(HydrationModelStatus *s) { hydrationModelStatus.reset(s); }

    void updateYourself(TimeStep *tStep);
    void printOutputAt(FILE *file, TimeStep *tStep);
};

class HydrationModelInterface : public Interface
{
    // interface for material access to HydrationModel
protected:
    std :: unique_ptr< HydrationModel > hydrationModel;
    double castAt;
    double constantHydrationDegree;

public:
    HydrationModel *giveHydrationModel() { return hydrationModel.get(); }
    virtual ~HydrationModelInterface() {}
    IRResultType initializeFrom(InputRecord *ir);

    contextIOResultType saveContext(DataStream &stream, ContextMode mode, void *obj = NULL)
    {
        if ( hydrationModel ) {
            hydrationModel->saveContext(stream, mode, obj);
        }

        return CIO_OK;
    }
    contextIOResultType restoreContext(DataStream &stream, ContextMode mode, void *obj = NULL)
    {
        if ( hydrationModel ) {
            hydrationModel->restoreContext(stream, mode, obj);
        }

        return CIO_OK;
    }

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