latticetransmat.C

Go to the documentation of this file.

/*
 *
 *                 #####    #####   ######  ######  ###   ###
 *               ##   ##  ##   ##  ##      ##      ## ### ##
 *              ##   ##  ##   ##  ####    ####    ##  #  ##
 *             ##   ##  ##   ##  ##      ##      ##     ##
 *            ##   ##  ##   ##  ##      ##      ##     ##
 *            #####    #####   ##      ######  ##     ##
 *
 *
 *             OOFEM : Object Oriented Finite Element Code
 *
 *               Copyright (C) 1993 - 2025   Borek Patzak
 *
 *
 *
 *       Czech Technical University, Faculty of Civil Engineering,
 *   Department of Structural Mechanics, 166 29 Prague, Czech Republic
 *
 *  This library is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU Lesser General Public
 *  License as published by the Free Software Foundation; either
 *  version 2.1 of the License, or (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *  Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public
 *  License along with this library; if not, write to the Free Software
 *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 */
 
#include "tm/Materials/LatticeMaterials/latticetransmat.h"
#include "domain.h"
#include "gausspoint.h"
#include "tm/Elements/LatticeElements/latticetransportelement.h"
#include "mathfem.h"
#include "staggeredproblem.h"
#include "classfactory.h"
#ifdef __SM_MODULE
 #include "sm/Elements/LatticeElements/latticestructuralelement.h"
#endif
 
namespace oofem {
REGISTER_Material(LatticeTransportMaterial);
 
void
LatticeTransportMaterial :: initializeFrom(InputRecord &ir)
{
    Material :: initializeFrom(ir);
 
    IR_GIVE_FIELD(ir, this->viscosity, _IFT_LatticeTransportMaterial_vis);
 
    IR_GIVE_FIELD(ir, this->permeability, _IFT_LatticeTransportMaterial_k);
 
    IR_GIVE_FIELD(ir, this->thetaS, _IFT_LatticeTransportMaterial_thetas);
 
    this->thetaR = 0.;
    IR_GIVE_OPTIONAL_FIELD(ir, this->thetaR, _IFT_LatticeTransportMaterial_thetar);
 
    //Options are 0 = constant conductivity and capacity and 1 = van Genuchten conductivity and capactity
    conType = 0;
    IR_GIVE_OPTIONAL_FIELD(ir, this->conType, _IFT_LatticeTransportMaterial_contype);
    this->capacity = 0;
    if ( conType == 0 ) {
        IR_GIVE_OPTIONAL_FIELD(ir, this->capacity, _IFT_LatticeTransportMaterial_c);
    } else if ( conType == 1 ) {
        IR_GIVE_FIELD(ir, this->paramM, _IFT_LatticeTransportMaterial_m);
 
        IR_GIVE_FIELD(ir, this->paramA, _IFT_LatticeTransportMaterial_a);
 
        //Assume that original van Genuchten is used.
        this->thetaM = this->thetaS;
        IR_GIVE_OPTIONAL_FIELD(ir, this->thetaM, _IFT_LatticeTransportMaterial_thetam);
 
        this->suctionAirEntry = 0;
        IR_GIVE_OPTIONAL_FIELD(ir, this->suctionAirEntry, _IFT_LatticeTransportMaterial_paev);
 
        if ( thetaM < thetaS ) {
            OOFEM_WARNING("thetaM cannot be smaller than thetaS. Choose thetaM=thetaS.");
            thetaM = thetaS;
        }
 
        //Define value of either air entry pressure or thetaM for later use. Only based on input parameters.
        if ( suctionAirEntry == 0. ) {
            suctionAirEntry = paramA * pow( ( pow( ( thetaS - thetaR ) / ( thetaM - thetaR ), -1. / paramM ) - 1. ), ( 1 - paramM ) );
        } else {
            thetaM = ( thetaS - thetaR ) * pow(1. + pow(suctionAirEntry / paramA, 1. / ( 1. - paramM ) ), paramM) + thetaR;
        }
    } else if ( conType != 0 && conType != 1 ) {
        OOFEM_ERROR("unknown conType mode");
    }
 
    crackTortuosity = 1.;
    IR_GIVE_OPTIONAL_FIELD(ir, this->crackTortuosity, _IFT_LatticeTransportMaterial_ctor);
 
    crackLimit = -1.;
    IR_GIVE_OPTIONAL_FIELD(ir, this->crackLimit, _IFT_LatticeTransportMaterial_clim);
}
 
 
FloatArrayF< 3 >
LatticeTransportMaterial :: computeFlux3D(const FloatArrayF< 3 > &grad, double field, GaussPoint *gp, TimeStep *tStep) const
{
    auto status = static_cast< LatticeTransportMaterialStatus * >( this->giveStatus(gp) );
    status->setTempGradient(grad);
    status->setTempField(field);
 
    double c = this->computeConductivity(field, gp, tStep);
    auto flux = -c * grad;
    status->setTempFlux(flux);
    return flux;
}
 
 
double
LatticeTransportMaterial :: give(int aProperty, GaussPoint *gp) const
{
    if ( aProperty == 'k' ) {
        return permeability;
    } else if ( ( aProperty == HeatCapaCoeff ) || ( aProperty == 'c' ) ) {
        return ( this->give('d', gp) );
    }
 
    return this->Material :: give(aProperty, gp);
}
 
 
double
LatticeTransportMaterial :: giveCharacteristicValue(MatResponseMode mode,
                                                    GaussPoint *gp,
                                                    TimeStep *tStep) const
{
    auto status = static_cast< LatticeTransportMaterialStatus * >( this->giveStatus(gp) );
    double suction = status->giveTempField();
 
    if ( mode == Capacity ) {
        return computeCapacity(suction, gp);
    } else if ( mode == Conductivity ) {
        return computeConductivity(suction, gp, tStep);
    } else {
        OOFEM_ERROR("unknown mode");
    }
 
    //return 0; // to make compiler happy
}
 
 
double
LatticeTransportMaterial :: computeConductivity(double suction,
                                                GaussPoint *gp,
                                                TimeStep *tStep) const
{
    auto status = static_cast< LatticeTransportMaterialStatus * >( this->giveStatus(gp) );
 
    double density = this->give('d', gp);
 
    double relativePermeability = 0.;
    double conductivity = 0.;
 
    double saturation, partOne, partTwo, numerator, denominator;
    if ( suction < this->suctionAirEntry || conType == 0 ) {
        saturation = 1.;
        relativePermeability = 1.;
    } else {
        partOne = pow(suction / this->paramA, 1. / ( 1. - this->paramM ) );
 
        saturation = ( this->thetaM - this->thetaR ) / ( this->thetaS - this->thetaR ) * pow(1. + partOne, -this->paramM);
 
        partTwo = ( this->thetaS - this->thetaR ) / ( this->thetaM - this->thetaR );
 
        numerator = ( 1. - pow(1. - pow(partTwo * saturation, 1 / this->paramM), this->paramM) );
        denominator = ( 1. - pow(1. - pow(partTwo, 1 / this->paramM), this->paramM) );
 
        relativePermeability = sqrt(saturation) * pow( ( numerator ) / ( denominator ), 2.0 );
    }
 
    //Calculate mass for postprocessing
    double mass = ( saturation * ( this->thetaS - this->thetaR ) + this->thetaR ) * density;
 
    status->setMass(mass);
 
    conductivity = this->permeability / this->viscosity * relativePermeability;
 
 
 
    //add crack contribution;
    //Read in crack lengths
    FloatArray crackLengths;
    static_cast< LatticeTransportElement * >( gp->giveElement() )->giveCrackLengths(crackLengths);
    FloatArray crackWidths;
    crackWidths.resize(crackLengths.giveSize() );
 
#ifdef __SM_MODULE
    IntArray coupledModels;
    if ( domain->giveEngngModel()->giveMasterEngngModel() ) {
        ( static_cast< StaggeredProblem * >( domain->giveEngngModel()->giveMasterEngngModel() ) )->giveCoupledModels(coupledModels);
        int couplingFlag = ( static_cast< LatticeTransportElement * >( gp->giveElement() ) )->giveCouplingFlag();
 
        if ( couplingFlag == 1 && coupledModels.at(1) != 0 && !tStep->isTheFirstStep() ) {
            IntArray couplingNumbers;
            static_cast< LatticeTransportElement * >( gp->giveElement() )->giveCouplingNumbers(couplingNumbers);
            for ( int i = 1; i <= crackLengths.giveSize(); i++ ) {
                if ( couplingNumbers.at(i) != 0 ) {
                    crackWidths.at(i) = static_cast< LatticeStructuralElement * >( domain->giveEngngModel()->giveMasterEngngModel()->giveSlaveProblem(coupledModels.at(1) )->giveDomain(1)->giveElement(couplingNumbers.at(i) ) )->giveCrackWidth();
                } else {
                    crackWidths.at(i) = 0.;
                }
            }
        }
    }
#endif
    //Read in crack widths from transport element
    if ( !domain->giveEngngModel()->giveMasterEngngModel() ) {
        static_cast< LatticeTransportElement * >( gp->giveElement() )->giveCrackWidths(crackWidths);
    }
    //Use crack width and apply cubic law
    double crackContribution = 0.;
 
    for ( int i = 1; i <= crackLengths.giveSize(); i++ ) {
        if ( crackWidths.at(i) < this->crackLimit || this->crackLimit < 0. ) {
            crackContribution += pow(crackWidths.at(i), 3.) / crackLengths.at(i);
        } else {
            printf("Limit is activated\n");
            crackContribution += pow(crackLimit, 3.) / crackLengths.at(i);
        }
    }
 
    crackContribution *=  this->crackTortuosity * relativePermeability / ( 12. * this->viscosity );
    conductivity += crackContribution;
    return density * conductivity;
}
 
 
double
LatticeTransportMaterial :: computeCapacity(double suction, GaussPoint *gp) const
{
    double cap = 0.;
 
    double density = this->give('d', gp);
 
    if ( conType == 0 ) {
        cap = this->capacity;
    } else {
        if ( suction < this->suctionAirEntry ) {
            cap = 0.;
        } else {
            double partOne = this->paramM / ( this->paramA * ( 1. - this->paramM ) );
            double partTwo = pow(suction / this->paramA, this->paramM / ( 1. - this->paramM ) );
            double partThree = pow(1. + pow(suction / this->paramA, 1. / ( 1. - this->paramM ) ), -this->paramM - 1.);
            cap = ( this->thetaM - this->thetaR ) * partOne * partTwo * partThree;
        }
    }
 
    return density * cap;
}
 
 
std::unique_ptr<MaterialStatus> 
LatticeTransportMaterial :: CreateStatus(GaussPoint *gp) const
{
    return std::make_unique<LatticeTransportMaterialStatus>(gp);
}
 
 
void
LatticeTransportMaterialStatus :: printOutputAt(FILE *File, TimeStep *tStep) const
{
    MaterialStatus :: printOutputAt(File, tStep);
 
    fprintf(File, "  state %.4e", field);
    fprintf(File, "  mass %.8e", mass);
    fprintf(File, "\n");
}
 
void
LatticeTransportMaterialStatus :: updateYourself(TimeStep *tStep)
{
    TransportMaterialStatus :: updateYourself(tStep);
}
 
 
void
LatticeTransportMaterialStatus :: initTempStatus()
{
    TransportMaterialStatus :: initTempStatus();
    oldPressure = field;
}
 
 
LatticeTransportMaterialStatus :: LatticeTransportMaterialStatus(GaussPoint *g) :
    TransportMaterialStatus(g)
{
    mass = 0.;
}
} // end namespace oofem