The isotropic damage models are based on the simplifying assumption
that the stiffness degradation is isotropic, i.e., stiffness moduli
corresponding to different directions decrease proportionally and
independently of direction of loading. Consequently, the damaged
stiffness matrix is expressed as

where is elastic stiffness matrix of the undamaged material and is the damage parameter. Initially, is set to zero, representing the virgin undamaged material, and the response is linear-elastic. As the material undergoes the deformation, the initiation and propagation of microdefects decreases the stiffness, which is represented by the growth of the damage parameter . For , the stiffness completely disappears.

In the present context, the
matrix represents the secant
stiffness that relates the total strain to the total stress

Similarly to the theory of plasticity, a loading function is introduced. In the damage theory, it is natural to work in the strain space and therefore the loading function is depending on the strain and on an additional parameter , describing the evolution of the damage. Physically, is a scalar measure of the largest strain level ever reached. The loading function usually has the form

where is the equivalent strain, i.e., the scalar measure of the strain level. Damage can grow only if current state reaches the boundary of elastic domain (). This is expressed by the following loading/unloading conditions

It remains to link the variable to the damage parameter . As both and grow monotonically, it is convenient to postulate an explicit evolution law

The important advantage of this explicit formulation is that the stress corresponding to the given strain can be evaluated directly, without the need to solve the nonlinear system of equations. For the given strain, the corresponding stress is computed simply by evaluating the current equivalent strain, updating the maximum previously reached equivalent strain value and the damage parameter and reducing the effective stress according to .

This general framework for computing stresses and
stiffness matrix is common for all material models of this type.
Therefore, it is natural to introduce
the base class for all isotropic-based damage models which provides the general
implementation for the stress and stiffness matrix evaluation
algorithms. The particular models then only provide their equivalent
strain and damage evolution law definitions.
The base class only declares the virtual services for computing equivalent
strain and corresponding damage. The implementation of common services
uses these virtual functions, but they are only declared at
**IsotropicDamageMaterial** class level and have to be
implemented by the derived classes.

Together with the material model, the corresponding status has to be
defined, containing all necessary history variables.
For the isotropic-based damage models, the only history variable is
the value of the largest strain level ever reached ().
In addition, the corresponding damage level will be stored.
This is not necessary because damage can be always computed from
corresponding .
The **IsotropicDamageMaterialStatus** class is derived from
**StructuralMaterialStatus** class. The base class represents the
base material status class for all structural statuses. At
**StructuralMaterialStatus** level, the attributes common to all
``structural analysis'' material models - the strain and
stress vectors (both the temporary and non-temporary) are introduced. The
corresponding services for accessing, setting, initializing, and
updating these attributes are provided.
Therefore, only the and parameters are introduced
(both the temporary and non-temporary). The corresponding services for
manipulating these attributes are added and services for context
initialization, update, and store/restore operations are overloaded, to
handle the history parameters properly.

class IsotropicDamageMaterialStatus : public StructuralMaterialStatus { protected: /// scalar measure of the largest strain level ever reached in material double kappa; /// non-equilibrated scalar measure of the largest strain level double tempKappa; /// damage level of material double damage; /// non-equilibrated damage level of material double tempDamage; public: /// Constructor IsotropicDamageMaterialStatus (int n, Domain*d, GaussPoint* g) ; /// Destructor ~IsotropicDamageMaterialStatus (); /// Prints the receiver state to stream void printOutputAt (FILE *file, TimeStep* tStep) ; /// Returns the last equilibrated scalar measure /// of the largest strain level double giveKappa () {return kappa;} /// Returns the temp. scalar measure of the /// largest strain level double giveTempKappa () {return tempKappa;} /// Sets the temp scalar measure of the largest /// strain level to given value void setTempKappa (double newKappa) { tempKappa = newKappa;} /// Returns the last equilibrated damage level double giveDamage () {return damage;} /// Returns the temp. damage level double giveTempDamage () {return tempDamage;} /// Sets the temp damage level to given value void setTempDamage (double newDamage) { tempDamage = newDamage;} // definition char* giveClassName (char* s) const { return strcpy(s,"IsotropicDamageMaterialModelStatus") ;} classType giveClassID () const { return IsotropicDamageMaterialStatusClass; } /** Initializes the temporary internal variables, describing the current state according to previously reached equilibrium internal variables. */ virtual void initTempStatus (); /** Update equilibrium history variables according to temp-variables. Invoked, after new equilibrium state has been reached. */ virtual void updateYourself(TimeStep*); // saves current context(state) into stream /** Stores context of receiver into given stream. Only non-temp internal history variables are stored. @param stream stream where to write data @param obj pointer to integration point, which invokes this method @return nonzero if o.k. */ contextIOResultType saveContext (FILE* stream, void *obj = NULL); /** Restores context of receiver from given stream. @param stream stream where to read data @param obj pointer to integration point, which invokes this method @return nonzero if o.k. */ contextIOResultType restoreContext(FILE* stream, void *obj = NULL); };

The base **IsotropicDamageMaterial** class is derived
from the **StructuralMaterial** class.

/** Base class representing general isotropic damage model. It is based on isotropic damage concept, assuming that damage evolution law is postulated in explicit form, relating damage parameter (omega) to scalar measure of the largest strain level ever reached in material (kappa). */ class IsotropicDamageMaterial : public StructuralMaterial { protected : /// Coefficient of thermal dilatation double tempDillatCoeff; /// Reference to bulk (undamaged) material LinearElasticMaterial *linearElasticMaterial; /** variable controling type of loading/unloading law, default set to idm_strainLevel defines the two possibilities: - idm_strainLevelCR the unloaing takes place, when strain level is smaller than the largest level ever reached; - idm_damageLevelCR the unloaing takes place, when damage level is smaller than the largest damage ever reached; */ enum loaUnloCriterium { idm_strainLevelCR, idm_damageLevelCR } llcriteria; public : /// Constructor IsotropicDamageMaterial (int n,Domain* d) ; /// Destructor ~IsotropicDamageMaterial () ; /// Returns nonzero indicating that receiver is nonlinear int hasNonLinearBehaviour () { return 1 ;} /** Tests, if material supports material mode. @param mode required material mode @return nonzero if supported, zero otherwise */ int hasMaterialModeCapability (MaterialMode mode); char* giveClassName (char* s) const { return strcpy(s,"IsotropicDamageMaterial") ;} classType giveClassID () const { return IsotropicDamageMaterialClass;} /// Returns reference to undamaged (bulk) material LinearElasticMaterial* giveLinearElasticMaterial () {return linearElasticMaterial;} /** Computes full 3d material stiffness matrix at given integration point, time, respecting load history in integration point. @param answer computed results @param form material response form @param mode material response mode @param gp integration point @param atTime time step (most models are able to respond only when atTime is current time step) */ virtual void give3dMaterialStiffnessMatrix (FloatMatrix& answer, MatResponseForm form, MatResponseMode mode, GaussPoint* gp, TimeStep* atTime); /** Computes the real stress vector for given strain increment and integration point. Temp vars are updated accordingly @param answer contains result @param form material response form @param gp integration point @param reducedStrain strain vector in reduced form @param tStep current time step (most models are able to respond only when tStep is current time step) */ void giveRealStressVector (FloatArray& answer, MatResponseForm form, GaussPoint* gp, const FloatArray &reducedStrain, TimeStep* tStep); /** Returns the integration point corresponding value in Reduced form. @param answer contain corresponding ip value, zero sized if not available @param aGaussPoint integration point @param type determines the type of internal variable @returns nonzero if ok, zero if var not supported */ virtual int giveIPValue (FloatArray& answer, GaussPoint* aGaussPoint, InternalStateType type, TimeStep* atTime) ; /** Returns the mask of reduced indexes of Internal Variable component . @param answer mask of Full VectorSize, with components being the indexes to reduced form vectors. @param type determines the internal variable requested @returns nonzero if ok or error is generated for unknown mat mode. */ virtual int giveIntVarCompFullIndx (IntArray& answer, InternalStateType type, MaterialMode mmode); /** Returns the type of internal variable (scalar, vector, ...) @param type determines the type of internal variable @returns type of internal variable */ virtual InternalStateValueType giveIPValueType (InternalStateType type); /** Returns the corresponding integration point value size in Reduced form. @param type determines the type of internal variable @returns var size, zero if var not supported */ virtual int giveIPValueSize (InternalStateType type, GaussPoint* aGaussPoint) ; /** Returns a vector of coefficients of thermal dilatation in direction of each material principal (local) axis. @param answer vector of thermal dilatation coefficients @param gp integration point @param tStep time step (most models are able to respond only when atTime is current time step) */ virtual void giveThermalDilatationVector (FloatArray& answer, GaussPoint*, TimeStep*) ; /** Computes the equivalent strain measure from given strain vector (full form). @param kappa return param, comtaining the corresponding equivalent strain @param strain total strain vector in full form @param gp integration point @param atTime timestep */ virtual void computeEquivalentStrain (double& kappa, const FloatArray& strain, GaussPoint* gp, TimeStep* atTime) = 0; /** computes the value of damage parameter omega, based on given value of equivalent strain @param omega contains result @param kappa equivalent strain measure */ virtual void computeDamageParam (double& omega, double kappa, const FloatArray& strain, GaussPoint* gp) =0; /** Instanciates the receiver from input record */ IRResultType initializeFrom (InputRecord* ir); protected: /** Abstract service allowing to perfom some initialization, when damage first appear @param kappa scalar measure of strain level @param totalStrainVector current total strain vector @param gp integration point */ virtual void initDamaged (double kappa, FloatArray& totalStrainVector, GaussPoint* gp) {} /// Creates corresponding material status MaterialStatus* CreateStatus (GaussPoint *gp) {return new IsotropicDamageMaterialStatus (1,domain,gp);} // Overloaded to use specialized versions of these // services possibly implemented by linearElastic member /** Method for computing plane stress stifness matrix of receiver. Default implementation overloaded to use direct implementation of corresponding service at bulk material model level. @param answer stifness matrix @param form material response form @param mode material response mode @param gp integration point, which load history is used @param atTime time step (most models are able to respond only when atTime is current time step) */ void givePlaneStressStiffMtrx (FloatMatrix& answer, MatResponseForm form, MatResponseMode, GaussPoint* gp, TimeStep* atTime); /** Method for computing plane strain stifness matrix of receiver. Default implementation overloaded to use direct implementation of corresponding service at bulk material model level. @param answer stifness matrix @param form material response form @param mode material response mode @param gp integration point, which load history is used @param atTime time step (most models are able to respond only when atTime is current time step) */ void givePlaneStrainStiffMtrx (FloatMatrix& answer, MatResponseForm form, MatResponseMode, GaussPoint* gp, TimeStep* atTime); /** Method for computing 1d stifness matrix of receiver. Default implementation overloaded to use direct implementation of corresponding service at bulk material model level. @param answer stifness matrix @param form material response form @param mode material response mode @param gp integration point, which load history is used @param atTime time step (most models are able to respond only when atTime is current time step) */ void give1dStressStiffMtrx (FloatMatrix& answer, MatResponseForm form, MatResponseMode, GaussPoint* gp, TimeStep* atTime); } ;

#include "isodamagemodel.h" #include "gausspnt.h" #include "flotmtrx.h" #include "flotarry.h" #include "structuralcrosssection.h" #include "mathfem.h" IsotropicDamageMaterial :: IsotropicDamageMaterial (int n, Domain *d) : StructuralMaterial (n,d) // // constructor // { linearElasticMaterial = NULL; llcriteria = idm_strainLevelCR; } IsotropicDamageMaterial :: ~IsotropicDamageMaterial () // // destructor // { delete linearElasticMaterial; } int IsotropicDamageMaterial :: hasMaterialModeCapability (MaterialMode mode) // // returns whether receiver supports given mode // { if ((mode == _3dMat) || (mode == _PlaneStress) || (mode == _PlaneStrain) || (mode == _1dMat)) return 1; return 0; }

Next, the implementation for core services is presented. They compute material stiffness matrices for various material modes and compute the stress vector. The later service computes the real stress corresponding to the previous history and the given strain vector. It computes a new locally consistent (equilibrated) state, which is stored in temporary variables of the corresponding status (attribute of the integration point).

void IsotropicDamageMaterial :: give3dMaterialStiffnessMatrix (FloatMatrix& answer, MatResponseForm form, MatResponseMode mode, GaussPoint* gp, TimeStep* atTime) // // computes full constitutive matrix for case of gp stress-strain state. // { IsotropicDamageMaterialStatus *status = (IsotropicDamageMaterialStatus*) this -> giveStatus (gp); double om = status->giveTempDamage(); om = min (om, 0.999999); this->giveLinearElasticMaterial()-> give3dMaterialStiffnessMatrix (answer, form, mode, gp, atTime); answer.times (1.0-om); return ; } void IsotropicDamageMaterial :: giveRealStressVector (FloatArray& answer, MatResponseForm form, GaussPoint* gp, const FloatArray& totalStrain, TimeStep* atTime) // // returns real stress vector in 3d stress space of receiver according to // previous load history and given strain. // After the new local point equilibrium is reached all the // temp history variables are updated to new local equilibrium state. // // { IsotropicDamageMaterialStatus *status = (IsotropicDamageMaterialStatus*) this -> giveStatus (gp); LinearElasticMaterial* lmat = this->giveLinearElasticMaterial(); FloatArray strainVector, reducedTotalStrainVector; FloatMatrix de; double f, equivStrain, tempKappa, omega; this->initGpForNewStep(gp); // substract stress independent part // note: eigenStrains (tepmerature) is not contained in mechanical // strain stored in gp therefore it is necessary to substract // always the total eigen strain value this->giveStressDependentPartOfStrainVector(reducedTotalStrainVector, gp, totalStrain,atTime, TotalMode); // compute equivalent strain this->computeEquivalentStrain (equivStrain, reducedTotalStrainVector, gp, atTime); if (llcriteria == idm_strainLevelCR) { // compute value of loading function if strainLevel crit apply f = equivStrain - status->giveKappa(); if (f <= 0.0) { // damage do not grow tempKappa = status->giveKappa(); omega = status->giveDamage(); } else { // damage grow tempKappa = equivStrain; this->initDamaged (tempKappa, reducedTotalStrainVector, gp); // evaluate damage parameter this->computeDamageParam (omega, tempKappa, reducedTotalStrainVector, gp); } } else if (llcriteria == idm_damageLevelCR) { // evaluate damage parameter first tempKappa = equivStrain; this->initDamaged (tempKappa, reducedTotalStrainVector, gp); this->computeDamageParam (omega, tempKappa, reducedTotalStrainVector, gp); if (omega < status->giveDamage()) { // unloading takes place omega = status->giveDamage(); } } else _error ("giveRealStressVector: unsupported loading/uloading criteria"); lmat->giveCharacteristicMatrix(de, ReducedForm, SecantStiffness, gp, atTime); de.times(1.0-omega); answer.beProductOf (de, reducedTotalStrainVector); // update gp status-> letTempStrainVectorBe (totalStrain); status-> letTempStressVectorBe (answer); status-> setTempKappa (tempKappa); status-> setTempDamage(omega); return ; } /* The following routines are not necessary, since the base StructuralMaterial class implements general algorithm how to derive the stiffness and compliance matrices for specific modes from full 3d stiffness/compliance. However, these general services are not efficient, since they require matrix inversion. The better and efficient way is to overload these service and provide efficient implementation, if available. */ void IsotropicDamageMaterial::givePlaneStressStiffMtrx (FloatMatrix& answer, MatResponseForm form, MatResponseMode mode, GaussPoint* gp, TimeStep* atTime) { IsotropicDamageMaterialStatus *status = (IsotropicDamageMaterialStatus*) this -> giveStatus (gp); double om = status->giveTempDamage(); om = min (om, 0.999999); this->giveLinearElasticMaterial()-> giveCharacteristicMatrix (answer, form, mode, gp, atTime); answer.times (1.0-om); return ; } void IsotropicDamageMaterial::givePlaneStrainStiffMtrx (FloatMatrix& answer, MatResponseForm form, MatResponseMode mode, GaussPoint* gp, TimeStep* atTime) { IsotropicDamageMaterialStatus *status = (IsotropicDamageMaterialStatus*) this -> giveStatus (gp); double om = status->giveTempDamage(); om = min (om, 0.999999); this->giveLinearElasticMaterial()-> giveCharacteristicMatrix (answer, form, mode, gp, atTime); answer.times (1.0-om); return ; } void IsotropicDamageMaterial::give1dStressStiffMtrx (FloatMatrix& answer, MatResponseForm form, MatResponseMode mode, GaussPoint* gp, TimeStep* atTime) { IsotropicDamageMaterialStatus *status = (IsotropicDamageMaterialStatus*) this -> giveStatus (gp); double om = status->giveTempDamage(); om = min (om, 0.999999); this->giveLinearElasticMaterial()-> giveCharacteristicMatrix (answer, form, mode, gp, atTime); answer.times (1.0-om); return ; }

int IsotropicDamageMaterial::giveIPValue (FloatArray& answer, GaussPoint* aGaussPoint, InternalStateType type, TimeStep* atTime) { IsotropicDamageMaterialStatus* status = (IsotropicDamageMaterialStatus*) this -> giveStatus (aGaussPoint); if (type == IST_DamageTensor) { answer.resize(1); answer.at(1) = status->giveDamage(); return 1; } else if (type == IST_DamageTensorTemp) { answer.resize(1); answer.at(1) = status->giveTempDamage(); return 1; } else if (type == IST_MaxEquivalentStrainLevel) { answer.resize(1); answer.at(1) = status->giveKappa (); return 1; }else return StructuralMaterial::giveIPValue (answer, aGaussPoint, type, atTime); } InternalStateValueType IsotropicDamageMaterial::giveIPValueType (InternalStateType type) { if ((type == IST_DamageTensor)) return ISVT_TENSOR; else return StructuralMaterial::giveIPValueType (type); } int IsotropicDamageMaterial::giveIntVarCompFullIndx (IntArray& answer, InternalStateType type, MaterialMode mmode) { if ((type == IST_DamageTensor)) { answer.resize (9); answer.at(1) = 1; return 1; } else return StructuralMaterial::giveIntVarCompFullIndx (answer, type, mmode); } int IsotropicDamageMaterial::giveIPValueSize (InternalStateType type, GaussPoint* aGaussPoint) { if ((type == IST_DamageTensor)) return 1; else return StructuralMaterial::giveIPValueSize (type, aGaussPoint); } void IsotropicDamageMaterial :: giveThermalDilatationVector (FloatArray& answer, GaussPoint * gp, TimeStep* tStep) // // returns a FloatArray(6) of initial strain vector // eps_0 = {exx_0, eyy_0, ezz_0, gyz_0, gxz_0, gxy_0}^T // caused by unit temperature in direction of // gp (element) local axes // { answer.resize (6); answer.zero(); answer.at(1) = this->tempDillatCoeff; answer.at(2) = this->tempDillatCoeff; answer.at(3) = this->tempDillatCoeff; return ; } IRResultType IsotropicDamageMaterial :: initializeFrom (InputRecord* ir) { const char *__proc = "initializeFrom"; // Required by IR_GIVE_FIELD macro IRResultType result; // Required by IR_GIVE_FIELD macro IR_GIVE_FIELD (ir, tempDillatCoeff, IFT_IsotropicDamageMaterial_talpha, "talpha"); // Macro return StructuralMaterial::initializeFrom (ir); }

In the last part, the implementation of
the **IsotropicDamageMaterialStatus** class services follows.
The *initTempStatus* and *updateYourself* services are
used to initialize temporary history variables according to the previous
equilibrium state or to update the non-temporary variables, when a new
equilibrium is reached. Also, the implementation of
*storeContext* and *restoreContext* services is shown.
The service *printYourself* is also overloaded. It allows to print history variables
into output file.
Note that all overloaded methods call the corresponding method of the base
class (to ensure that processing corresponding to the base class is
performed, to guarantee consistency) and then the required functionality is added.

IsotropicDamageMaterialStatus::IsotropicDamageMaterialStatus (int n, Domain* d, GaussPoint* g) :StructuralMaterialStatus(n,d,g) { kappa = tempKappa = 0.0; damage = tempDamage = 0.0; } IsotropicDamageMaterialStatus::~IsotropicDamageMaterialStatus () {} void IsotropicDamageMaterialStatus :: printOutputAt (FILE *file, TimeStep* tStep) { StructuralMaterialStatus :: printOutputAt (file, tStep); fprintf (file,"status { "); if (this->damage > 0.0) { fprintf (file,"kappa %f, damage %f ",this->kappa, this->damage); } fprintf (file,"}\n"); } void IsotropicDamageMaterialStatus::initTempStatus () { StructuralMaterialStatus :: initTempStatus(); this->tempKappa = this->kappa; this->tempDamage= this->damage; } void IsotropicDamageMaterialStatus::updateYourself(TimeStep* atTime) { StructuralMaterialStatus::updateYourself(atTime); this->kappa = this->tempKappa; this->damage= this->tempDamage; } contextIOResultType IsotropicDamageMaterialStatus::saveContext (FILE* stream, void *obj) { contextIOResultType iores; // save parent class status if ((iores = StructuralMaterialStatus :: saveContext (stream, obj)) != CIO_OK) THROW_CIOERR(iores); // write a raw data if (fwrite(&kappa,sizeof(double),1,stream) != 1) THROW_CIOERR(CIO_IOERR); if (fwrite(&damage,sizeof(double),1,stream)!= 1) THROW_CIOERR(CIO_IOERR); return CIO_OK; } contextIOResultType IsotropicDamageMaterialStatus::restoreContext(FILE* stream, void *obj) { contextIOResultType iores; // read parent class status if ((iores = StructuralMaterialStatus :: restoreContext(stream,obj)) != CIO_OK) THROW_CIOERR(iores); // read raw data if (fread (&kappa,sizeof(double),1,stream) != 1) THROW_CIOERR(CIO_IOERR); if (fread (&damage,sizeof(double),1,stream) != 1) THROW_CIOERR(CIO_IOERR); return CIO_OK; }

The above general class provides the general framework. The derived
classes have to implement only services
*computeEquivalentStrain* and *computeDamageParam*,
which are only related to particular material model under consideration.

Borek Patzak 2018-01-02