Material models¶
The base class for all material models is the Material class, derived from the FEMComponent parent. It declares analysis independent part of the material interface. The analysis specific part of the interface should be introduced by derived classes, which are supposed to represent the base classes for specific problem. The typical example is StructuralMaterial class which declares the services of material model necessary for the structural analysis. The generic services declared or implemented at top Material level include
Material status related services. The material model has to be able to create instance of corresponding material status, where the history variables are stored. This is done by invoking CreateStatus. The status corresponding to a given integration point can be requested using the giveStatus service.
Services for integration point update and initialization. There are generally two sets of history variables kept in corresponding material statuses for each integration point. One set is referring to previous equilibrated state, the second one to the actual state during the solution (more precisely to the achieved local-equilibrium state), which may not correspond to the global equilibrium state. The methods are provided to update the actual state as equilibrated (updateYourself service) and for initialization of actual state to previous equilibrium (initTempStatus). The implementation of these services simply extract the corresponding status of a given integration point and calls the corresponding service of material status.
Services for testing the material model capabilities (testMaterialExtension, hasMaterialModeCapability, and hasNonLinearBehaviour services).
Services for requesting internal variables and properties.
Analysis specific material model classes¶
The direct derived classes of the Material class, such as StructuralMaterial are supposed to declare the analysis specific part of the material model interface, which is required (and assumed) by the corresponding element class (StructuralElement) and cross section class (StructuralCrossSection).
To store all necessary history variables of the material model, so called status concept is adopted. The material status can be thought as container of all necessary history variables. Usually two kinds of these variables are stored. The temporary ones refer to the actual state of an integration point, but do not necessary correspond to the global equilibrium. These are changing during global equilibrium search iteration. The non-temporary variables are related to the previously converged state. For each material model, the corresponding twin status has to be defined, and the unique instance for each integration point has to be created and associated with it. The integration point provides the services for accessing corresponding status. All material statuses, related to particular material models, have to be derived from the base MaterialStatus class. This class declares the basic status interface. The two most important services are:
initTempStatus intended to initialize the temporary internal variables according to variables related to previously reached equilibrium state,
updateYourself designed to update the equilibrium-like history variables according to temporary variables, when the new global equilibrium has been reached. The derived classes should also define methods for accessing the corresponding history variables.
Structural Material class - Example¶
The structural (or mechanical) constitutive model should generally support several so-called material models, corresponding to various modeling assumptions (plane-stress, plane-strain, or 1D-behavior, for example). The concept of multiple material modes is generally supported by Material class, which provides the services for testing the model capabilities. It is generally assumed, that results obtained from constitutive model services are formed according to material mode. This mode is attribute of each integration point, which is compulsory parameter of all material services. For computational convenience, the so-called full and reduced formats of stress/strains vectors are introduced, corresponding to material modes. The full format includes all components, even if they are zero due to stress/strain mode nature. In the reduced format, only generally nonzero components are stored. (Full format should be used only if absolutely necessary, to avoid wasting of space. For example, it is used by output routines to print results in general form). Methods for converting vectors between full and reduced format are provided. If possible, all computations should be performed in reduced space.
Important
The convention used to construct reduced strain/stress vectors is as follows: If in a particular mode a particular stress component is zero, the corresponding strain is not computed and not stored in reduced vector, and in full vector there is zero value on corresponding position. On the other hand, if zero strain component is imposed, then this condition must be taken into account in geometrical relations (at element level), and corresponding components are included in stress/strain reduced vectors.
Generally, the following major tasks are declared by StructuralMaterial or inherited from Material class:
Computing the real stress vector (tensor) at an integration point for a given strain increment and updating its temporary state corresponding to the local equilibrium, but not necessarily to the global equilibrium (see giveRealStressVector). The parameters include the total strain vector and the corresponding integration point. The total strain is defined as strain computed directly from the displacement field at a given time. The stress independent parts (temperature, eigen strains) should be subtracted and the corresponding load-history variables (stored in the corresponding status) can be used. The temporary history variables in the status should be updated according to the newly reached state. The temporary history variables are moved into equilibrium history variables just after the global structure equilibrium has been reached by the iteration process.
Updating the integration point state (final state), when the global equilibrium has been reached.
Returning material stiffness and/or flexibility matrices for a given material mode. The general methods computing the response for the specific material mode are provided, based on converting 3D stiffness or compliance matrix to that corresponding to the specific material mode. But, if it is possible to compute stiffness or compliance matrix directly for the specific mode, then these general methods should be overloaded.
Storing/restoring integration point state to/from a stream.
Requesting internal variables, their types and properties (giveIPValue, giveIPValueSize, giveIPValueType, and giveIntVarCompFullIndx services).
Returning material properties.
Structural material services should not be called directly by elements. Instead, the elements should pass their requests to the corresponding cross section model, that performs all necessary integration over its volume and invokes corresponding material model services.
The StructuralMaterial class comes with definition of associated material status - StructuralMaterialStatus. This is only an abstract class. For every instance of StructuralMaterial class there should be a twin statu class, which defines and maintains all necessary state variables. The StructuralMaterialStatus only introduces attributes common to all structural analysis material models - the strain and stress vectors (referring both to temporary state, corresponding to the local equilibrium and non-temporary ones, corresponding to the global equilibrium). The corresponding services for accessing, setting, initializing, and updating these attributes are provided.
Isotropic Damage Model - Example¶
In this section, the implementation of ann isotropic damage model will be described. To cover the various models based on isotropic damage concept, a base class IsotropicDamageMaterial is defined first, declaring the necessary services and providing the implementation of them, which are general. The derived classes then only implement a particular damage-evolution law.
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 \({D}_e\) is elastic stiffness matrix of the undamaged material and \(\omega\) is the damage parameter. Initially, \(\omega\) 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 \(\omega\). For \(\omega=1\), the stiffness completely disappears.
In the present context, the \(D\) matrix represents the secant stiffness that relates the total strain to the total stress
Similarly to the theory of plasticity, a loading function $f 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 \(\kappa\), describing the evolution of the damage. Physically, \(\kappa\) is a scalar measure of the largest strain level ever reached. The loading function usually has the form
where \(\tilde\varepsilon\) 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 (\(f=0\)). This is expressed by the following loading/unloading conditions
It remains to link the variable \(\kappa\) to the damage parameter \(\omega\). As both \(\kappa\) and \(\omega\) 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 \(\kappa\) and the damage parameter and reducing the effective stress according to \({\sigma} = (1-\omega){D}_e {\varepsilon}\).
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 (\(\kappa\)). In addition, the corresponding damage level $omega$ will be stored. This is not necessary because damage can be always computed from corresponding \(\kappa\). 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 \(\kappa\) and \(\omega\) 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.
The declaration of the IsotropicDamageMaterialStatus class follows:
1 class IsotropicDamageMaterialStatus : public StructuralMaterialStatus {
2 protected:
3 /// scalar measure of the largest strain level ever reached in material
4 double kappa;
5 /// non-equilibrated scalar measure of the largest strain level
6 double tempKappa;
7 /// damage level of material
8 double damage;
9 /// non-equilibrated damage level of material
10 double tempDamage;
11
12 public:
13 /// Constructor
14 IsotropicDamageMaterialStatus (int n, Domain*d, GaussPoint* g) ;
15 /// Destructor
16 ~IsotropicDamageMaterialStatus ();
17
18 /// Prints the receiver state to stream
19 void printOutputAt (FILE *file, TimeStep* tStep) ;
20
21 /// Returns the last equilibrated scalar measure
22 /// of the largest strain level
23 double giveKappa () {return kappa;}
24 /// Returns the temp. scalar measure of the
25 /// largest strain level
26 double giveTempKappa () {return tempKappa;}
27 /// Sets the temp scalar measure of the largest
28 /// strain level to given value
29 void setTempKappa (double newKappa) { tempKappa = newKappa;}
30 /// Returns the last equilibrated damage level
31 double giveDamage () {return damage;}
32 /// Returns the temp. damage level
33 double giveTempDamage () {return tempDamage;}
34 /// Sets the temp damage level to given value
35 void setTempDamage (double newDamage) { tempDamage = newDamage;}
36
37
38 // definition
39 char* giveClassName (char* s) const
40 { return strcpy(s,"IsotropicDamageMaterialModelStatus") ;}
41 classType giveClassID () const
42 { return IsotropicDamageMaterialStatusClass; }
43
44 /**
45 Initializes the temporary internal variables,
46 describing the current state according to
47 previously reached equilibrium internal variables.
48 */
49 virtual void initTempStatus ();
50 /**
51 Update equilibrium history variables
52 according to temp-variables.
53 Invoked, after new equilibrium state has been reached.
54 */
55 virtual void updateYourself(TimeStep*);
56
57 // saves current context(state) into stream
58 /**
59 Stores context of receiver into given stream.
60 Only non-temp internal history variables are stored.
61 @param stream stream where to write data
62 @param obj pointer to integration point, which invokes this method
63 @return nonzero if o.k.
64 */
65 contextIOResultType saveContext (FILE* stream, void *obj = NULL);
66 /**
67 Restores context of receiver from given stream.
68 @param stream stream where to read data
69 @param obj pointer to integration point, which invokes this method
70 @return nonzero if o.k.
71 */
72 contextIOResultType restoreContext(FILE* stream, void *obj = NULL);
73 };
The base IsotropicDamageMaterial class is derived from the StructuralMaterial class. The generic methods defined by StructuralMaterial require only to implement two fundamental services:
give3dMaterialStiffnessMatrix for evaluating the material stiffness matrix in full, 3D mode,
giveRealStressVector_3d for computing the real stress vector in 3D mode.
The remaining material modes are derived from these two basic services. However, the derived classes can reimplement (overload) these services, when desired. In our case, it is possible to obtain reduced stress and stiffness in more efficient way, than done in generic way on StructuralMaterial level and which requires matrix inversion. Therefore, the services for reduced modes are overloaded.
1 class IsotropicDamageMaterial : public StructuralMaterial
2 {
3 protected:
4 /// Coefficient of thermal dilatation.
5 double tempDillatCoeff = 0.;
6
7 /// Maximum limit on omega. The purpose is elimination of a too compliant material which may cause convergence problems. Set to something like 0.99 if needed.
8 double maxOmega = 0.999999;
9
10 /// Indicator of the type of permanent strain formulation (0 = standard damage with no permanent strain)
11 int permStrain = 0;
12
13 /// Reference to bulk (undamaged) material
14 LinearElasticMaterial *linearElasticMaterial = nullptr;
15 /**
16 * Variable controlling type of loading/unloading law, default set to idm_strainLevel
17 * defines the two two possibilities:
18 * - idm_strainLevelCR the unloading takes place, when strain level is smaller than the largest level ever reached;
19 * - idm_damageLevelCR the unloading takes place, when damage level is smaller than the largest damage ever reached;
20 */
21 enum loaUnloCriterium { idm_strainLevelCR, idm_damageLevelCR } llcriteria = idm_strainLevelCR;
22
23 public:
24 /// Constructor
25 IsotropicDamageMaterial(int n, Domain *d);
26 /// Destructor
27 virtual ~IsotropicDamageMaterial();
28
29 bool hasMaterialModeCapability(MaterialMode mode) const override;
30 const char *giveClassName() const override { return "IsotropicDamageMaterial"; }
31
32 /// Returns reference to undamaged (bulk) material
33 LinearElasticMaterial *giveLinearElasticMaterial() { return linearElasticMaterial; }
34
35 FloatMatrixF<6,6> give3dMaterialStiffnessMatrix(MatResponseMode mode, GaussPoint *gp, TimeStep *tStep) const override;
36
37 void giveRealStressVector(FloatArray &answer, GaussPoint *gp,
38 const FloatArray &reducedStrain, TimeStep *tStep) override;
39
40 FloatArrayF<6> giveRealStressVector_3d(const FloatArrayF<6> &strain, GaussPoint *gp, TimeStep *tStep) const override
41 {
42 FloatArray answer;
43 const_cast<IsotropicDamageMaterial*>(this)->giveRealStressVector(answer, gp, strain, tStep);
44 return answer;
45 }
46 FloatArrayF<4> giveRealStressVector_PlaneStrain( const FloatArrayF<4> &strain, GaussPoint *gp, TimeStep *tStep) const override
47 {
48 FloatArray answer;
49 const_cast<IsotropicDamageMaterial*>(this)->giveRealStressVector(answer, gp, strain, tStep);
50 return answer;
51 }
52 FloatArray giveRealStressVector_StressControl(const FloatArray &strain, const IntArray &strainControl, GaussPoint *gp, TimeStep *tStep) const override
53 {
54 FloatArray answer;
55 const_cast<IsotropicDamageMaterial*>(this)->giveRealStressVector(answer, gp, strain, tStep);
56 return answer;
57 }
58 FloatArrayF<3> giveRealStressVector_PlaneStress(const FloatArrayF<3> &strain, GaussPoint *gp, TimeStep *tStep) const override
59 {
60 FloatArray answer;
61 const_cast<IsotropicDamageMaterial*>(this)->giveRealStressVector(answer, gp, strain, tStep);
62 return answer;
63 }
64 FloatArrayF<1> giveRealStressVector_1d(const FloatArrayF<1> &strain, GaussPoint *gp, TimeStep *tStep) const override
65 {
66 FloatArray answer;
67 const_cast<IsotropicDamageMaterial*>(this)->giveRealStressVector(answer, gp, strain, tStep);
68 return answer;
69 }
70
71 int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep) override;
72
73 FloatArrayF<6> giveThermalDilatationVector(GaussPoint *gp, TimeStep *tStep) const override;
74 virtual double evaluatePermanentStrain(double kappa, double omega) const { return 0.; }
75
76 /**
77 * Returns the value of material property 'aProperty'. Property must be identified
78 * by unique int id. Integration point also passed to allow for materials with spatially
79 * varying properties
80 * @param aProperty ID of property requested.
81 * @param gp Integration point,
82 * @return Property value.
83 */
84 double give(int aProperty, GaussPoint *gp) const override;
85 /**
86 * Computes the equivalent strain measure from given strain vector (full form).
87 * @param[out] kappa Return parameter, containing the corresponding equivalent strain.
88 * @param strain Total strain vector in full form.
89 * @param gp Integration point.
90 * @param tStep Time step.
91 */
92 virtual double computeEquivalentStrain(const FloatArray &strain, GaussPoint *gp, TimeStep *tStep) const = 0;
93 /**Computes derivative of the equivalent strain with regards to strain
94 * @param[out] answer Contains the resulting derivative.
95 * @param strain Strain vector.
96 * @param gp Integration point.
97 * @param tStep Time step.
98 */
99 virtual void computeEta(FloatArray &answer, const FloatArray &strain, GaussPoint *gp, TimeStep *tStep) const { OOFEM_ERROR("not implemented"); }
100 /**
101 * Computes the value of damage parameter omega, based on given value of equivalent strain.
102 * @param[out] omega Contains result.
103 * @param kappa Equivalent strain measure.
104 * @param strain Total strain in full form.
105 * @param gp Integration point.
106 */
107 virtual double computeDamageParam(double kappa, const FloatArray &strain, GaussPoint *gp) const = 0;
108
109 void initializeFrom(InputRecord &ir) override;
110 void giveInputRecord(DynamicInputRecord &input) override;
111
112 MaterialStatus *CreateStatus(GaussPoint *gp) const override { return new IsotropicDamageMaterialStatus(gp); }
113
114 FloatMatrixF<1,1> give1dStressStiffMtrx(MatResponseMode mmode, GaussPoint *gp,
115 TimeStep *tStep) const override;
116
117 void saveContext(DataStream &stream, ContextMode mode) override;
118 void restoreContext(DataStream &stream, ContextMode mode) override;
119
120 protected:
121 /**
122 * Abstract service allowing to perform some initialization, when damage first appear.
123 * @param kappa Scalar measure of strain level.
124 * @param totalStrainVector Current total strain vector.
125 * @param gp Integration point.
126 */
127 virtual void initDamaged(double kappa, FloatArray &totalStrainVector, GaussPoint *gp) const { }
128
129 /**
130 * Returns the value of derivative of damage function
131 * wrt damage-driving variable kappa corresponding
132 * to a given value of the kappa, depending on
133 * the type of selected damage law.
134 * @param kappa Equivalent strain measure.
135 * @param gp Integration point.
136 */
137 virtual double damageFunctionPrime(double kappa, GaussPoint *gp) const {
138 OOFEM_ERROR("not implemented");
139 return 0;
140 }
141
142 FloatMatrixF<3,3> givePlaneStressStiffMtrx(MatResponseMode mmode, GaussPoint *gp,TimeStep *tStep) const override;
143 FloatMatrixF<4,4> givePlaneStrainStiffMtrx(MatResponseMode mmode, GaussPoint *gp, TimeStep *tStep) const override;
144
145 };
Finally we present the implementation of the selected methods of IsotropicDamageMaterial class. Let us start with the implementation of the give3dMaterialStiffnessMatrix method, which computes the 3D material stiffness, respecting the previous loading history described by state variable:
1 FloatMatrixF<6,6>
2 IsotropicDamageMaterial :: give3dMaterialStiffnessMatrix(MatResponseMode mode,
3 GaussPoint *gp,
4 TimeStep *tStep) const
5 //
6 // computes full constitutive matrix for case of gp stress-strain state.
7 //
8 {
9 auto status = static_cast< IsotropicDamageMaterialStatus * >( this->giveStatus(gp) );
10 double tempDamage;
11 if ( mode == ElasticStiffness ) {
12 tempDamage = 0.0;
13 } else {
14 tempDamage = status->giveTempDamage();
15 tempDamage = min(tempDamage, maxOmega);
16 }
17
18 auto d = this->linearElasticMaterial->give3dMaterialStiffnessMatrix(mode, gp, tStep);
19 return d * (1.0 - tempDamage);
20 //TODO - correction for tangent mode
21 }
The giveRealStressesVector method computes the real stress respecting the previous loading history (described by state variables) 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).
1 void
2 IsotropicDamageMaterial :: giveRealStressVector(FloatArray &answer, GaussPoint *gp,
3 const FloatArray &totalStrain,
4 TimeStep *tStep)
5 //
6 // returns real stress vector in 3d stress space of receiver according to
7 // previous level of stress and current
8 // strain increment, the only way, how to correctly update gp records
9 //
10 {
11 IsotropicDamageMaterialStatus *status = static_cast< IsotropicDamageMaterialStatus * >( this->giveStatus(gp) );
12 LinearElasticMaterial *lmat = this->giveLinearElasticMaterial();
13 FloatArray reducedTotalStrainVector;
14 FloatMatrix de;
15 double f, equivStrain, tempKappa = 0.0, omega = 0.0;
16
17 this->initTempStatus(gp);
18
19 // subtract stress-independent part
20 // note: eigenStrains (temperature) are present in strains stored in gp
21 // therefore it is necessary to subtract always the total eigen strain value
22 this->giveStressDependentPartOfStrainVector(reducedTotalStrainVector, gp, totalStrain, tStep, VM_Total);
23
24 // compute equivalent strain
25 equivStrain = this->computeEquivalentStrain(reducedTotalStrainVector, gp, tStep);
26
27 if ( llcriteria == idm_strainLevelCR ) {
28 // compute value of loading function if strainLevel crit apply
29 f = equivStrain - status->giveKappa();
30
31 if ( f <= 0.0 ) {
32 // damage does not grow
33 tempKappa = status->giveKappa();
34 omega = status->giveDamage();
35 } else {
36 // damage grows
37 tempKappa = equivStrain;
38 this->initDamaged(tempKappa, reducedTotalStrainVector, gp);
39 // evaluate damage parameter
40 omega = this->computeDamageParam(tempKappa, reducedTotalStrainVector, gp);
41 }
42 } else if ( llcriteria == idm_damageLevelCR ) {
43 // evaluate damage parameter first
44 tempKappa = equivStrain;
45 this->initDamaged(tempKappa, reducedTotalStrainVector, gp);
46 omega = this->computeDamageParam(tempKappa, reducedTotalStrainVector, gp);
47 if ( omega < status->giveDamage() ) {
48 // unloading takes place
49 omega = status->giveDamage();
50 }
51 } else {
52 OOFEM_ERROR("unsupported loading/unloading criterion");
53 }
54
55 // get material stiffness from bulk material
56 lmat->giveStiffnessMatrix(de, SecantStiffness, gp, tStep);
57
58 // damage deactivation in compression for 1D model
59 if ( ( reducedTotalStrainVector.giveSize() > 1 ) || ( reducedTotalStrainVector.at(1) > 0. ) ) {
60 //emj
61 de.times(1.0 - omega);
62 }
63
64 answer.beProductOf(de, reducedTotalStrainVector);
65
66 // update gp status
67 status->letTempStrainVectorBe(totalStrain);
68 status->letTempStressVectorBe(answer);
69 status->setTempKappa(tempKappa);
70 status->setTempDamage(omega);
71 }
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. For full reference, please refer to the src/sm/Materials/isotropicdamagematerial.h and src/sm/Materials/isotropicdamagematerial.C files.
For example of simple Isotropic Damage Model, please refer to the src/sm/Materials/ConcreteMaterials/idm1.h and src/sm/Materials/ConcreteMaterials/idm1.C files.