How to use the Constitutive Law class

From KratosWiki
(Difference between revisions)
Jump to: navigation, search
(Conventions)
 
(28 intermediate revisions by 3 users not shown)
Line 14: Line 14:
  
 
''voigt notation:''
 
''voigt notation:''
 +
 
- 3D case:
 
- 3D case:
 
   '''STRAIN''' Voigt Notation:  e00 e11 e22 2*e01 2*e12 2*e02
 
   '''STRAIN''' Voigt Notation:  e00 e11 e22 2*e01 2*e12 2*e02
Line 22: Line 23:
 
   '''STRESS''' Voigt Notation:  s00 s11 s22  s01   
 
   '''STRESS''' Voigt Notation:  s00 s11 s22  s01   
  
- 2D plane stress
+
- 2D plane stress/strain
 
   '''STRAIN''' Voigt Notation:  e00 e11 2*e01  
 
   '''STRAIN''' Voigt Notation:  e00 e11 2*e01  
 
   '''STRESS''' Voigt Notation:  s00 s11  s01
 
   '''STRESS''' Voigt Notation:  s00 s11  s01
  
 +
The constitutive law works on the basis of the '''total deformation gradient F''', defined as
 +
  F  := D(X) / D(X0)
 +
that is, as the deformation gradient connecting the original and deformed configuration
  
The constitutive law API is also designed to ease the implementation not only of Total Lagrangian but also of Updated Lagrandian and "Spatial Lagrangian" approaches.
+
where the initial position X0 is the one obtained by
In order to understant the interface we shall consider 3 possible references:
+
  const array_1d<double,3>& X0 = node->GetInitialPosition()
 +
and the deformed one by
 +
  const array_1d<double,3>& X = node->Coordinates()
 +
  //must coincide with      X = node->GetInitialPosition() + node.FastGetSolutionStepValue(DISPLACEMENT);
  
1 - The "Initial Reference", representing the position in space occupied by the nodes at the very beginning of the simulation, before any deformation is applied. We recall that in Kratos such position can be recovered by
+
The ConstitutiveLaw '''always returns the total stress'''. Formulations expressed in terms of strain increments shall store internally the strain stresses from which the increment
    const array_1d<double,3>& X0 = node->InitialCoordinates()
+
shall be computed
2 - The "last known position", identifying the position in space occupied in the last moment at which everything is considered known. We could name such position as "Xlast"
+
3 - The "final" configuration of the mesh, identified as the position at which the nodes are located at the very latest iteration. This position can be obtained as
+
        const array_1d<double,3>& X = node->Coordinates()
+
        //must coincide with          node->InitialCoordinates() + node.FastGetSolutionStepValue(DISPLACEMENT);
+
  
taking into account this definitions, the Constitutive Law defines
+
== Usage API ==
 +
The constitutive law API is based on the use of an auxiliary "Parameters" data structure, designed to encapsulate the data to be passed to the CL and received from it.
 +
The parameters data structure should be initialized using the following constructor:
 +
      Parameters (const GeometryType& rElementGeometry
 +
                  ,const Properties& rMaterialProperties
 +
                  ,const ProcessInfo& rCurrentProcessInfo)
  
- F0 := D(Xref) / D(X0)
+
Thus allowing to encapsulate the pointer to the elemental properties, to the element geometry and to the process info.
+
- F  := D(X) / D(Xref)
+
  
The total deformation gradient can hence be always obtained as
+
The data structure '''does not contain any internal storage''' and should be initialized with pointers to memory ''owned by the caller element''.
  Ftot = prod(F,F0)
+
Full documentation of the code can be found in the file [https://kratos.cimne.upc.es/projects/kratos/repository/entry/kratos/kratos/includes/constitutive_law.h constitutive_law.h].
 +
For ease, the getter interface, ''returning a reference to the encapsulated data'', is reported here
  
taking into account such definitions the CL API allows the user to define on which domain configuration to compute the given stress/strain measure
+
      GetOptions() //returns a reference to a flag container, to be employed in passing options to the CL
There are thus 3 different possibilities
+
      GetDeterminantF() 
 +
      GetDeformationGradientF()
 +
      GetShapeFunctionsValues()
 +
      GetShapeFunctionsDerivatives()
 +
      GetStrainVector() //INPUT/OUTPUT -- note that F will be used preferentially instead of the input strain
 +
      GetStressVector()
 +
      GetConstitutiveMatrix()
 +
      GetProcessInfo()
 +
      GetMaterialProperties()
 +
      GetElementGeometry()
  
1 - INITIAL CONFIGURATION - in this case everything should be computed on the reference configuration, so that
+
The "Options" flag represents the fundamental tool in steering the control of the constitutive law behaviour. The interface provides a number of boolean flags that can be passed to the constitutive law:
     F0 := I
+
 
     F  := D(X) / D(X0)  
+
    COMPUTE_STRAIN
     In this case the user shall instruct the constitutive law to obtain this behaviour by setting
+
    COMPUTE_STRESS
 +
    COMPUTE_CONSTITUTIVE_TENSOR
 +
 
 +
ask the constitutive law to compute these magnitudes and will be returned in the variables of the Data structure GetStrainVector(), GetStressVector() and GetConstitutiveMatrix() respectively.
 +
 
 +
    COMPUTE_STRAIN_ENERGY
 +
 
 +
asks the constitutive law to compute the strain energy, elastic or plastic, after that you can get it using a GetValue(kratos<variable>, double).
 +
 
 +
    ISOCHORIC_TENSOR_ONLY     
 +
    VOLUMETRIC_TENSOR_ONLY   
 +
    TOTAL_TENSOR               
 +
    FINALIZE_MATERIAL_RESPONSE
 +
 
 +
these flags are used internally in the constitutive law to set the type of calculation needed depending on the law and the particular time sequence.
 +
 
 +
 
 +
A fundamental feature of the constitutive law is to implement internally the transformations between different stress measures.
 +
When a user writes an element he should decide what
 +
stress measure is desired. The list of available options is provided in the enum:
 +
    enum StressMeasure
 +
    {
 +
        StressMeasure_PK1,            //stress related to reference configuration non-symmetric
 +
        StressMeasure_PK2,           //stress related to reference configuration
 +
        StressMeasure_Kirchhoff,      //stress related to current  configuration
 +
        StressMeasure_Cauchy          //stress related to current  configuration
 +
     }
 +
 
 +
at the moment of writing the element the developer should hence query to the constitutive law for the desired stress measure. This is achieved by picking
 +
one of the functions:
 +
     CalculateMaterialResponsePK1( parameters )
 +
    CalculateMaterialResponsePK2( parameters )
 +
     CalculateMaterialResponseKirchhoff( parameters )
 +
    CalculateMaterialResponseCauchy( parameters )
 +
 
 +
the elasticity tensor and the corresponding stresses will be stored in the "parameters" and shall be accessed by the Getter funtions described above.
 +
 
 +
At the end of a given solution step (typically in the FinalizeSolutionStep function of the element), internal variables should be updated by calling the function
 +
    FinalizeMaterialResponsePK1( parameters )
 +
    FinalizeMaterialResponsePK2( parameters )
 +
    FinalizeMaterialResponseKirchhoff( parameters )
 +
    FinalizeMaterialResponseCauchy( parameters )
 +
 
 +
An example of usage of the ConstitutiveLaw from within a total lagrangian element could be as follows:
 +
    Parameters parameters(GetGeometry(), GetProperties(), rCurrentProcessInfo);
 
      
 
      
2 - LAST KNOWN CONFIGURATION - in this case everything should be computed on the reference configuration, so that
+
    //here we essentially set the input parameters
     F0 := D(Xref) / D(X0) ---> RICCARDOS PROPOSAL: should change this to Fref
+
    parameters.SetDeterminantF(detF) //assuming the determinant is computed somewhere else
     F  := D(X) / D(Xref) ---> RICCARDOS PROPOSAL: should change the name to DF
+
    parameters.SetDeformationGradientF(F) //F computed somewhere else
3 - INITIAL CONFIGURATION - in this case everything should be computed on the reference configuration, so that
+
   
     F0 := D(X) / D(X0)  
+
    //here we set the space on which the results shall be written
    F := I
+
    Matrix ConstitutiveTensor(strain_size,strain_size); //note that C is allocated in the element
 +
     Vector stress(strain_size);
 +
    parameters.SetConstitutiveMatrix(ConstitutiveTensor) //assuming the determinant is computed somewhere else
 +
     parameters.SetStressVector(stress) //F computed somewhere else
 +
   
 +
    //instruct the constitutive law that both stress and the constitutivetensor are needed
 +
    parameters.GetOptions().Set(COMPUTE_STRESS, True)
 +
    parameters.GetOptions().Set(COMPUTE_CONSTITUTIVE_TENSOR, True)
 +
   
 +
    //actually do the computations in the ConstitutiveLaw   
 +
    constitutivelaw->CalculateMaterialResponsePK2(parameters); //here the calculations are actually done
 +
   
 +
    //here stress and C are already updated since we passed them to the CL
 +
 
 +
The constitutive law also provides helper functions to allow pushing/pulling stresses from one configuration to another.
 +
 
 +
== Providing the Deformation gradient F or the strain ==
 +
Within Kratos the strain can be either provided to the constitutive law as an input parameter, OR obtained internally within the constitutive law
 +
depending on the kinematics chosen.
 +
 
 +
=== OPTION 1 -- PROVIDE STRAIN AS INPUT PARAMETER ===
 +
This option is chosen by selecting
 +
    flags.Set(COMPUTE_STRAIN, false)
 +
 
 +
in this case the constitutive law reads the value of "StrainVector" and uses it as a basis for the computations of the stresses
 +
 
 +
=== OPTION 2 -- PROVIDE F AS INPUT PARAMETER ===
 +
The idea in this case is to provide a DeformationGradient '''F''' as an input parameter.
 +
 
 +
In the case of small deformation kinematics, is is not customary to define the deformation gradient. For this case (small deformations) one one can construct the deformation gradient given the strain vector as
 +
 
 +
  F(0,0) = 1+strain(0);     F(0,1) = 0.5*strain(3); F(0,2) = 0.5*strain(5)
 +
  F(1,0) = 0.5*strain(3);  F(1,1) = 1+strain(1);  F(1,2) = 0.5*strain(4)
 +
  F(2,0) = 0.5*strain(5);  F(2,1) = 0.5*strain(4); F(2,2) = 1+strain(2)
 +
 
 +
for the 3D case.
 +
 
 +
== Additional Material ==
 +
[[How to use 3D laws for PlaneStress cases]]
 +
 
 +
[[How to use constitutive law to describe RVE behaviour]]
 +
 
 +
[[How to use an Abaqus "umat" function through the interface]]

Latest revision as of 00:13, 19 March 2016

The constitutive law behaviour is dealt with in kratos by the use of the class "ConstitutiveLaw", with a public interface defined in the file

  kratos/kratos/includes/constitutive_law.h

which also provides some rather extensive inline documentation (in the form of comments in the code).

By design such file aims to provide a very flexible interface to constitutive law modelling, with the specific goal of maximizing the flexibility in the implementation of complex constitutive behaviours. While such approach provide obvious advantages, it also implies that the API is more complex than what would be strictly needed for very simple constitutive laws.

The objective of current HowTo is to provide a brief introduction to the interface

Contents

Conventions

Through the whole section, the following convenctions will be employed:

voigt notation:

- 3D case:

 STRAIN Voigt Notation:  e00 e11 e22 2*e01 2*e12 2*e02
 STRESS Voigt Notation:  s00 s11 s22   s01   s12   s02
        

- 2D plane strain/axisymmetric case (4 stress components)

 STRAIN Voigt Notation:  e00 e11 e22 2*e01 
 STRESS Voigt Notation:  s00 s11 s22   s01   

- 2D plane stress/strain

 STRAIN Voigt Notation:  e00 e11 2*e01 
 STRESS Voigt Notation:  s00 s11   s01

The constitutive law works on the basis of the total deformation gradient F, defined as

  F  := D(X) / D(X0) 

that is, as the deformation gradient connecting the original and deformed configuration

where the initial position X0 is the one obtained by

  const array_1d<double,3>& X0 = node->GetInitialPosition()

and the deformed one by

  const array_1d<double,3>& X = node->Coordinates() 
  //must coincide with      X = node->GetInitialPosition() + node.FastGetSolutionStepValue(DISPLACEMENT);

The ConstitutiveLaw always returns the total stress. Formulations expressed in terms of strain increments shall store internally the strain stresses from which the increment shall be computed

Usage API

The constitutive law API is based on the use of an auxiliary "Parameters" data structure, designed to encapsulate the data to be passed to the CL and received from it. The parameters data structure should be initialized using the following constructor:

     Parameters (const GeometryType& rElementGeometry
                 ,const Properties& rMaterialProperties
                 ,const ProcessInfo& rCurrentProcessInfo)

Thus allowing to encapsulate the pointer to the elemental properties, to the element geometry and to the process info.

The data structure does not contain any internal storage and should be initialized with pointers to memory owned by the caller element. Full documentation of the code can be found in the file constitutive_law.h. For ease, the getter interface, returning a reference to the encapsulated data, is reported here

     GetOptions() //returns a reference to a flag container, to be employed in passing options to the CL
     GetDeterminantF()   
     GetDeformationGradientF() 
     GetShapeFunctionsValues() 
     GetShapeFunctionsDerivatives() 
     GetStrainVector() //INPUT/OUTPUT -- note that F will be used preferentially instead of the input strain
     GetStressVector() 
     GetConstitutiveMatrix() 
     GetProcessInfo() 
     GetMaterialProperties() 
     GetElementGeometry()

The "Options" flag represents the fundamental tool in steering the control of the constitutive law behaviour. The interface provides a number of boolean flags that can be passed to the constitutive law:

   COMPUTE_STRAIN
   COMPUTE_STRESS
   COMPUTE_CONSTITUTIVE_TENSOR

ask the constitutive law to compute these magnitudes and will be returned in the variables of the Data structure GetStrainVector(), GetStressVector() and GetConstitutiveMatrix() respectively.

   COMPUTE_STRAIN_ENERGY

asks the constitutive law to compute the strain energy, elastic or plastic, after that you can get it using a GetValue(kratos<variable>, double).

   ISOCHORIC_TENSOR_ONLY      
   VOLUMETRIC_TENSOR_ONLY     
   TOTAL_TENSOR                
   FINALIZE_MATERIAL_RESPONSE

these flags are used internally in the constitutive law to set the type of calculation needed depending on the law and the particular time sequence.


A fundamental feature of the constitutive law is to implement internally the transformations between different stress measures. When a user writes an element he should decide what stress measure is desired. The list of available options is provided in the enum:

   enum StressMeasure
   {
       StressMeasure_PK1,            //stress related to reference configuration non-symmetric
       StressMeasure_PK2,            //stress related to reference configuration
       StressMeasure_Kirchhoff,      //stress related to current   configuration
       StressMeasure_Cauchy          //stress related to current   configuration
   }

at the moment of writing the element the developer should hence query to the constitutive law for the desired stress measure. This is achieved by picking one of the functions:

   CalculateMaterialResponsePK1( parameters )
   CalculateMaterialResponsePK2( parameters )
   CalculateMaterialResponseKirchhoff( parameters )
   CalculateMaterialResponseCauchy( parameters )

the elasticity tensor and the corresponding stresses will be stored in the "parameters" and shall be accessed by the Getter funtions described above.

At the end of a given solution step (typically in the FinalizeSolutionStep function of the element), internal variables should be updated by calling the function

   FinalizeMaterialResponsePK1( parameters )
   FinalizeMaterialResponsePK2( parameters )
   FinalizeMaterialResponseKirchhoff( parameters )
   FinalizeMaterialResponseCauchy( parameters )

An example of usage of the ConstitutiveLaw from within a total lagrangian element could be as follows:

   Parameters parameters(GetGeometry(), GetProperties(), rCurrentProcessInfo);
   
   //here we essentially set the input parameters
   parameters.SetDeterminantF(detF) //assuming the determinant is computed somewhere else
   parameters.SetDeformationGradientF(F) //F computed somewhere else
   
   //here we set the space on which the results shall be written
   Matrix ConstitutiveTensor(strain_size,strain_size); //note that C is allocated in the element
   Vector stress(strain_size);
   parameters.SetConstitutiveMatrix(ConstitutiveTensor) //assuming the determinant is computed somewhere else
   parameters.SetStressVector(stress) //F computed somewhere else
   
   //instruct the constitutive law that both stress and the constitutivetensor are needed
   parameters.GetOptions().Set(COMPUTE_STRESS, True)
   parameters.GetOptions().Set(COMPUTE_CONSTITUTIVE_TENSOR, True)
   
   //actually do the computations in the ConstitutiveLaw    
   constitutivelaw->CalculateMaterialResponsePK2(parameters); //here the calculations are actually done 
   
   //here stress and C are already updated since we passed them to the CL

The constitutive law also provides helper functions to allow pushing/pulling stresses from one configuration to another.

Providing the Deformation gradient F or the strain

Within Kratos the strain can be either provided to the constitutive law as an input parameter, OR obtained internally within the constitutive law depending on the kinematics chosen.

OPTION 1 -- PROVIDE STRAIN AS INPUT PARAMETER

This option is chosen by selecting

   flags.Set(COMPUTE_STRAIN, false)

in this case the constitutive law reads the value of "StrainVector" and uses it as a basis for the computations of the stresses

OPTION 2 -- PROVIDE F AS INPUT PARAMETER

The idea in this case is to provide a DeformationGradient F as an input parameter.

In the case of small deformation kinematics, is is not customary to define the deformation gradient. For this case (small deformations) one one can construct the deformation gradient given the strain vector as

  F(0,0) = 1+strain(0);     F(0,1) = 0.5*strain(3); F(0,2) = 0.5*strain(5)
  F(1,0) = 0.5*strain(3);   F(1,1) = 1+strain(1);   F(1,2) = 0.5*strain(4)
  F(2,0) = 0.5*strain(5);   F(2,1) = 0.5*strain(4); F(2,2) = 1+strain(2)

for the 3D case.

Additional Material

How to use 3D laws for PlaneStress cases

How to use constitutive law to describe RVE behaviour

How to use an Abaqus "umat" function through the interface

Personal tools
Categories