Table of contents

×Contents

User Material Models

Introduction

IMPETUS Solver comes with support for user defined material models. Our solution is to let the user compile a DLL (dynamic-link library) file which then is linked to the solver during the start of a simulation. The user material is not compiled into the software itself and it allows for several advantages. Since object code of the software isn't needed, it will works with any version of the solver (unless we've done major changes to the interface). This means one can use the latest version of the Solver. We allows up to five different user defined material models.

Structure

The user material zip file (provided by our support) contains a Visual Studio 2022 solution with two project folders inside. The C++ project "mat_user_functions" contains a wrapper to call built-in subroutines/functions from the solver. Do not edit the files inside this project. The other project "mat_user" which is an Intel Fortran project, contains the code for user defined material models. There are several files inside this project: mat_user.F, a few mat_XXXX.F files (examples) and custom_functions.F. As the name suggests, custom_functions.F can be used to define custom subroutines/functions. This is completely optional but we recommend using it. The file mat_user.F is the main file and this is the one that should be modified. Start by opening mat_user.sln with Visual Studio. Before compilation, remember to change compilation mode from Debug to Release.

Now, expand "mat_user" project from the solution tree and open mat_user.F under "Source Files". We've provided with some example models which are implemented in different files, but called by mat_user_1(), mat_user_2(), mat_user_3() and mat_user_4(). Feel free to modify these and customise them to your needs. The example in mat_user_1() is similar to our built-in *MAT_METAL material model.

After the compilation, mat_user.dll will be copied to the output folder located in the root directory of the Visual Studio project. When executing a simulation, check the "Use defined material model" checkbox and add the path to the file.

Defined subroutines in mat_user.F
mat_user_X()
Main subroutine for the user material X (X: 1-5)
init_mat_user_X()
Initialisation subroutine for the user material X (X: 1-5)
mat_user_set_state_variables()
Set the number of state variables, state variable output, texture location (for fibre directions), damage type location and temperature location.
Arguments for mat_user_X():
Argument
Description
strain(6)
Total strain tensor
dstrain(6)
Strain increment tensor
stress(6)
Stress tensor
F(9)
Deformation gradient
U(6)
Right stretch tensor
ifail
integer variable to set if the current integration point has failed or not
iface_update
Integer variable to force update of faces
hist(*)
History state variable array (Size defined in mat_user_set_state_variables())
cmat(100)
Material parameters
stiffness(4)
Used for time step calculation. 1 = shear, 2 = bulk, 3 = xi, 4 = bfac
iel
Element ID
ip
Integration point ID (local ID for the specific element)
itype
Element type
dt1
Current time step size
tt
Current time
Arguments for init_mat_user_X():
hist(*)
History state variable array. Size defined in mat_user_set_state_variables()
cmat(100)
Material parameters
stress(6)
Stress tensor
strain(6)
Strain tensor
stiffness(4)
Used for time step calculation. 1 = shear, 2 = bulk, 3 = xi, 4 = bfac
x
Node coordinate
iel
Element ID
ip
Integration point ID (local ID for the specific element)
itype
Element type

Getting started

Symmetric tensors

All symmetric tensors in our solver uses an alternative notation of the Voigt notation. The traditional voigt notation for a 3X3 matrix is defined as: 11, 22, 33, 23, 13, 12. Our definition is: 11, 22, 33, 12, 23, 13.

\begin{bmatrix} \sigma_{1} & \sigma_{4} & \sigma_{6} \\ & \sigma_{2} & \sigma_{5} \\ & & \sigma_{3} \end{bmatrix} \begin{bmatrix} \sigma_{1} & \sigma_{6} & \sigma_{5} \\ & \sigma_{2} & \sigma_{4} \\ & & \sigma_{3} \end{bmatrix}
Top matrix is our notation, bottom is the traditional notation.
Material parameters

Material parameters defined in the *MAT_USER_X command can be accessed by the cmat array. Position 1-6 are tied to the first row after the material ID definition (1st value). For example cmat(1) tied to the 2nd value, cmat(2) the 3rd, cmat(3) the 4th, etc. For a normal use case, we use the following: cmat(1) = density, cmat(2) = young's modulus, cmat(3) = poisson's ratio, etc.

Second row starts at index cmat(7), third row at cmat(15), forth row at cmat(23), etc.

There are a few slots in cmat which are reserved. cmat(43) and cmat(44) are reserved for pre-calculated shear and bulk modulus.

Example: material data from input: density, young's modulus and poisson's ratio
! *MAT_USER_1
! 1, 7800.0, 2.0e9, 0.3

dens  = cmat(1) ! 7800.0
young = cmat(2) ! 2.0e9
pr    = cmat(3) ! 0.3

! pre-calculated
shear = cmat(43)
bulk  = cmat(44)
Damage properties

If a damage command ID is defined in the input, then the damage properties will be stored in cmat(60) → cmat(70) and cmat(80) → cmat(93). 

! these variables are automatically set depending damage and erosion type
damage_type = nint(cmat(80))
erode_flg   = nint(cmat(81))

[...]

! damage evolution (see definition in the section "Built-in subroutines/functions")
if (damage_type.ne.0.and.(dmg.lt.1.0d0.or.erode_flg.ge.2)) then
  call damage(damage_type, erode_flg, cmat(60), cmat(82), dmg, dmg2,
              stress, sig_eff, p, strain, epsp, depsp, rate, T, evol, 0, 0, 0)
endif
Thermal properties

If a thermal command ID is defined in the input, then the thermal properties will be stored in cmat(71) → cmat(75). 

! thermal properties (see *PROP_THERMAL in manual for explanation)
hexp  = cmat(71)
Cp    = cmat(72)
hcond = cmat(73)
k     = cmat(74)
Tref  = cmat(75)

[...]

! update temperature T if we have plastic flow
if (depsp.gt.0.0d0.and.Cp.ne.0.0d0) then
  T = T + k*sig_eff*depsp/(Cp*dens)
endif
State variables

The state variable array is a container where one can store data for each integration point. These variables can change throughout the simulation. Examples of typical variables are "Effective Plastic Strain", "Damage" and temperature. To use state variables, use the hist array to load and store data. It's a good idea to set initial values for these variables in init_mat_user_X() subroutine.

Example of loading and storing state variables
! load state variables: effective plastic strain, damage, yield stress and temperature
epsp  = hist(1)
dmg   = hist(2)
sigy0 = hist(3)
T     = hist(4)

! some code involving epsp, dmg, sigy0 and/or T
[...]

! save state variables
hist(1) = epsp
hist(2) = dmg
hist(3) = sigy0
hist(4) = T

Before state variables can be used, one needs to define the size of the array. This is done by editing mat_user_set_state_variables() in mat_user.F. In our examples, these are set in the subsequential subroutine calls since they are implemented in different files). For example, see mat_template_state_variables() in mat_template.F (which is called by mat_user_set_state_variables() in mat_user.F).

For instance, if we are going to specify the number of history variables for mat_user_1(), we need to assign a value to num_variable(1). In the example (see code below), we have set the size to 3 for position 1 in the num_variable() array. Position 1 indicates user material ID 1, value 3 is the number of state variables for this user material.

Example
num_variable(1) = 4
State variables as contour plot

With the size set, we are ready to go but we will not be able to plot these values in the GUI. None of these values are written to our output files by default. To add a state variable to the output, we need to tell the solver to output it. This can be done by calling mat_user_set_name(). variable_id tells which index in the state variable array to output. variable_name specifies the name the state variable should have.

Example
user_material_id         = 1
variable_id              = 1
variable_name(1:40)      = "Effective plastic strain"
call mat_user_set_name(user_material_id, variable_id, variable_name)

variable_id              = 2
variable_name(1:40)      = "Damage"
call mat_user_set_name(user_material_id, variable_id, variable_name)

They can be found under the "General" category in Contour plot. The prefix "User: " will automatically be added to the name when using Contour plot in the GUI.

Curves & functions

Curves & functions can be used in the user material models just like in *MAT_METAL. To use curves, one needs to tell the solver where in the material parameter list the curve is located. This can be done by editing mat_user_set_state_variables() and modify/add the array user_curve_location. This is a 2D array with two arguments defining user curve ID and material ID.

Example of defining two curves
user_curve_location(1, 1) = 7
user_curve_location(2, 1) = 8

This will tell the solver that there are curves in positions 7 and 8 in the material parameter array (cmat) for mat_user_1. Position 7 and 8 represents the two first columns on the second line in the *MAT_USER_X command. More information about how to call curves and functions can be found in the section "Built-in subroutines/functions" further down in this document.

Texture location

When implementing material models using fibre direction (for example fibre composites) with *INITIAL_MATERIAL_DIRECTION_WRAP, a texture location must be set. The location tells IMPETUS Solver where in the state variables array (hist) to set the initial fibre direction. It requires 9 slots (3 vectors) in the state variables array (from the specified index) so make sure to allocate enough space when setting num_variable().

Example below shows allocation of 14 state variables and that we set the last 9 of them for initial fibre direction:

Example of allocation of 14 state variables, set the last 9 of them for initial fibre direction
num_variable(1)     = 14
texture_location(1) = 6
Temperature location

Temperature definition from for example *INITIAL_TEMPERATURE can be obtained through the temperature_location array. The value set in this array will tells IMPETUS Solver where in the state variables array (hist) to set the initial temperature.

Example of allocation of 14 state variables, set the 4th position for temperature
num_variable(1)         = 14
temperature_location(1) = 4
Damage definition

Material failure can be defined as element erosion or node splitting. The variable damage_location defines where in the state variable array that damage is located. This variable is only needed for node splitting but we recommend to set it anyway even when not using node splitting. erode_flg defines if the material model should use element erosion (=1), node splitting (=2) or none (=0). If erode_flg_pos is set, the value of the material parameter (cmat) in that specified location will determine to use element erosion or node splitting.

Example
damage_location          = 2
erode_flg                = 2
erode_flg_pos            = 0
call mat_user_set_damage(user_material_id, damage_location, erode_flg, erode_flg_pos)

This should enough for node splitting since the Solver will handle the rest. For element erosion, failure handling must be defined manually. To do this, add this after the damage calculations in mat_user_X():

Example
! element erosion (node splitting will not execute this if-statement)
if (dmg.ge.1.0d0.and.erode_flg.lt.2) then
dmg = 1.0d0

! reset deviatoric stresses
stress(:) = 0.0d0

! erosion -> set failure flag
if (erode_flg.eq.1) ifail = ifail + 1
endif

We set the dmg variable to 1.0, reset the stress array and we increment ifail with 1. If enough integration points per element fails, then the solver will automatically erode it. This code can co-exist with node splitting capability since node splitting will set the erode_flg variable to 2. If that happens then this code will never be executed.

Built-in subroutines/functions

We've made some of the core functionality of the solver available. Return values from the subroutines are marked with underscore in the parameter list. Functions return its value as return values.

Damage interface for *PROP_DAMAGE_X
call damage( type , erode_flg , dmg_fcn , dmg_parm , dmg , dmg2 , sig_dev , sig_eff , p , eps , epsp , depsp , rate , T , evol , Q , ctip , ip_vol )
Interface to invoke one of the built-in damage failure criterion models. The first four parameters should be passed on from cmat.
Parameter
Description
type
Damage type parameter (integer)
erode_flg
Erosion flag parameter (integer)
dmg_fcn
Use cmat(60) as parameter (double precision)
dmg_parm
Use cmat(82) as parameter (double precision)
dmg
Ductile damage parameter (double precision)
dmg2
Brittle damage parameter (double precision)
sig_dev
6-component deviatoric stress tensor (double precision)
sig_eff
Effective stress (double precision)
p
Pressure (double precision)
eps
6-component strain tensor (double precision)
epsp
Effective plastic strain (double precision)
depsp
Effective plastic strain increment (double precision)
rate
Effective strain rate (double precision)
T
Temperature (double precision)
evol
Element volume (double precision)
Q
3x3-matrix (element size and direction), used for regularisation (double precision)
ctip
Crack tip (experimental) (integer)
ip_vol
Integration point volume (double precision)
Example
damage_type = nint(cmat(80))
erode_flg   = nint(cmat(81))

[...]

! damage evolution
if (damage_type.ne.0.and.(dmg.lt.1.0d0.or.erode_flg.ge.2)) then
    call damage(damage_type, erode_flg, cmat(60), cmat(82), dmg, dmg2, stress,
                sig_eff, p, strain, epsp, depsp, rate, T, evol, 0, 0, 0)
endif
Load Curve / Function
call load_curve( idlc , x , f , p , epsp , sigy0 , T , dmg )
Returns the ordinate value of x in the defined idlc to the variable f. The last 5 parameters (p, epsp, sigy0, T & dmg) are for functions. They must be set to 0 if not being used. This call can be used both to call curves and functions.
See the command *CURVE and *FUNCTION for more information.
Parameter
Description
idlc
Function/Curve ID (integer)
x
Value of x in function f(x) (double precision)
f
Return value (double precision)
p
Pressure parameter (double precision)
epsp
Old effective plastic strain parameter (double precision)
sigy0
Yield stress parameter (double precision)
T
Thermal parameter (double precision)
dmg
Damage parameter (double precision)
Example
! Load a curve
call load_curve(idlc, epsp+deps, sigy1, 0, 0, 0, 0, 0)

! Load a curve or a function with pressure and effective strain rate variables set
call load_curve(idlc, epsp+deps, sigy1, p, epsp, 0, 0, 0)

! Load a curve or a function width all parameters set
call load_curve(idlc, epsp+deps, sigy1, p, epsp, sigy0, T, dmg)
Normalise a vector
call normalize( v , x )
This subroutine will return normalised v and also its length x.
Parameter
Description
v
3-component vector (double precision)
x
length (double precision)
Cross product
call cross( a1 , a2 , a3 )
This subroutine will return the cross product between a1 & a2 in a3.
Parameter
Description
a1
3-component vector (double precision)
a2
3-component vector (double precision)
a3
Return value, 3-component vector (double precision)
Dot product
value = dot( v1 , v2 )
Returns the dot product between v1 and v2.
Parameter
Description
v1
3-component vector (double precision)
v2
3-component vector (double precision)
Polar decomposition
call polar_decomposition( F , R , U , Uinv )
This subroutine computes the polar decomposition and returns new values for R, U & Uinv.
Parameter
Description
F
Positive definite tensor (3X3) (double precision)
R
Rotation tensor (3X3) (double precision)
U
Symmetric stretch tensor U (double precision)
Uinv
Inverted symmetric stretch tensor (double precision)
Eigenvalues of symmetric 3X3-matrix
call eigval_3x3_symm( avec , eigval )
This subroutine returns the eigen values calculated from avec.
Parameter
Description
avec
6-component tensor (double precision)
eigval
Returns 3-component vector (double precision)
Eigenvectors of symmetric 3X3-matrix
call eigvec_3x3_symm( avec , eigval , eigvec )
This subroutine returns the eigen vectors calculated from avec.
Parameter
Description
avec
6-component tensor (double precision)
eigval
3-component vector (double precision)
eigvec
Returns two eigenvectors eigvec(1) and eigvec(4) in one 6-component array (double precision)
Example
! First, find eigenvalues and then use them to find the eigenvectors
call eigval_3x3_symm(stress, eigval)
call eigvec_3x3_symm(stress, eigval, eigvec)
Invert 3X3-matrix
call inv_3x3( A , B , detA )
This subroutine reads matrix A and returns the inverted matrix to B. The determinate is also returned to detA.
Parameter
Description
A
3X3 matrix to be inverted (double precision)
B
Return value, inverted matrix (double precision)
detA
Return value, determinant (double precision)
Volume of an element
value = element_volume( element_id , element_type )
Returns volume of an element in double precision.
Parameter
Description
element_id
Element ID (integer)
element_type
Element type (integer)
Volume of an integration point
value = ip_volume( element_id , element_type , ip_id )
Returns volume of an integration point in double precision.
Parameter
Description
element_id
Element ID (integer)
element_type
Element type (integer)
ip_id
Integration point ID (integer)
Initial volume of an integration point
value = ip_initial_volume( element_id , element_type , ip_id )
Returns initial volume of an integration point in double precision.
Parameter
Description
element_id
Element ID (integer)
element_type
Element type (integer)
ip_id
Integration point ID (integer)
Get element nodes
call element_nodes( element_id , element_type , node_list )
This subroutine returns an array of node IDs of a specific element.
Parameter
Description
element_id
Element ID (integer)
element_type
Element type (integer)
node_list
Return value, array (integer)
Get element neighbours
call element_neighbor_list( element_id , element_type , element_list , num_neigh )
This subroutine returns an 2D array of neighbouring elements and the number of neighbours. First value of the array contains the neighbouring element ID. The second value contains the element type of the neighbouring element.
Parameter
Description
element_id
Element ID (integer)
element_type
Element type (integer)
element_list
Return value, 2D array (integer)
num_neigh
Return value, number of neighbours (integer)
Get integration point coordinate
call ip_coord( element_id , element_type , ip_id , ip_coord )
This subroutine returns the coordinate of an integration point.
Parameter
Description
element_id
Element ID (integer)
element_type
Element type (integer)
ip_id
Integration point ID (integer)
ip_coord
Return value, 3-component vector (double precision)
Element conversion (internal to user ID)
value = element_i2u( element_id , element_type )
Returns user (external) ID of an element.
Parameter
Description
element_id
Element ID (integer)
element_type
Element type (integer)