Using External Materials

Using External Materials

If you need to use a constitutive model which is not available among the built-in material models, it is possible to program it yourself. Such a material function, termed an external material, is coded in C. If you already have an existing code in another language like Fortran or C++, it is however possible to write a wrapper function to it.

Before moving to implementing your own material model, there are however two other options to consider:

•

Many material models like hyperelasticity, creep and plasticity have as one of the options in addition to the standard models. Any material model which you can describe using built-in variables is most conveniently described here.

•

A material model which can partially described in terms of a PDE can often be implemented using one of the mathematical interfaces. Stresses or strains computed in that interface are then injected in an existing material model using the External Stress and External Strain subnodes.

There are two basic types of external material functions: those which completely replace other material definitions in a domain, and those that just compute an inelastic strain contribution to be used as part of an existing material model. The former is referenced from an External Stress-Strain Relation node, whereas the latter is referenced from an External Strain subnode.

During the solution, an external material routine is always called for each Gauss point during evaluation of stiffness matrices and computation of residuals. During result presentation, the external material can be called from any location in the geometry, as requested by for example graphs and point evaluations.

Almost invariably, you need to store state variables in the external material, such as for example plastic strains. The state variables are stored at the Gauss points. If an external material is called at another location, the state variables will be interpolated to that location. This means that the state of the material may not be exactly consistent there, which can lead to some artifacts during result presentation. You can avoid this problem by using the gpeval operator.

In the COMSOL Multiphysics Reference Manual:

•

Working with External Materials

•

External Material

Library of Utility Functions

In order to simplify the task of writing the code for an external material, a library of utility routines is provided. It provides a toolkit for operations common in solid mechanics such as various tensor operations or computing principal stresses and equivalent stresses.

List of Utility Functions

csext_add: Function that adds two matrices

* Function: csext_add

* -------------------

* Adds two (3 x 3) matrices 'A' and 'B' and stores the

* result 'A + B' in 'C'.

* double A[3][3]

* double B[3][3]

* double C[3][3] (output)

* Return value:

void csext_add(double A[3][3],double B[3][3], double C[3][3]);

csext_det: Function that computes the matrix determinant

* Function: csext_det

* -------------------

* Computes the determinant of a (3 x 3) matrix 'A' and

* returns the value.

* double A[3][3]

* Return value:

double csext_det(double A[]3[3]);

csext_dev: Function that computes the deviator of a matrix

* Function: csext_dev

* -------------------

* Computes the deviator of a (3 x 3) matrix 'A'

* and stores the value in 'dev'.

* double A[3][3]

* double dev[3][3] (output)

* Return value:

void csext_dev(double A[3][3], double dev[3][3]);

csext_dot: Function that computes the inner product of two matrices

* Function: csext_dot

* -------------------

* Computes the inner product of (3 x 3) matrices 'A' and 'B',

* and returns 'A : B'.

* double A[3][3]

* double B[3][3]

* Return value:

double csext_dot(double A[3][3],double B[3][3]);

csext_eig: Function that computes the spectral decomposition of a symmetric matrix

* Function: csext_eig

* -------------------

* Computes the eigenvalues and eigenvectors of a

* symmetric (3 x 3) matrix 'A'.

* The eigenvalues are stored in 'vals', sorted with the

* largest value in vals[0].

* The eigenvectors are stored column-wise in 'vecs' with

* the same ordering as the eigenvalues.

* Normal execution returns 0. nonzero means that the

* computation failed.

* double A[3][3]

* double vals[3] (output)

* double vecs[3][3] (output)

* Return value:

int csext_eig(double A[3][3], double vals[3], double vecs[3][3]);

csext_gl: Function that computes the Green-Lagrange strain tensor

* Function: csext_gl

* ------------------

* Computes the Green-Lagrange strain tensor 'egl' based on

* the Right Cauchy-Green deformation tensor 'rcg'.

* double rcg[3][3]

* double egl[3][3] (output)

* Return value:

void csext_gl(double rcg[3][3],double egl[3][3]);

csext_inv: Function that computes the matrix inverse

* Function: csext_inv

* -------------------

* Computes the inverse of a (3 x 3) matrix 'A'.

* The inverse, if it exists, is stored in 'inv'.

* 0 is returned if successful, -1 if 'A' is numerically

* singular. The matrix 'A' is considered singular if

* abs(det(A))<tol.

* double A[3][3]

* double inv[3][3] (output)

* Return value:

int csext_inv(double A[3][3], double tol, double inv[3][3]);

csext_lcg: Function that computes the Left Cauchy–Green deformation tensor

* Function: csext_lcg

* -------------------

* Computes the Left Cauchy-Green deformation 'lcg' tensor

* based on the deformation gradient 'defgrad'.

* double defgrad[3][3]

* double lcg[3][3] (output)

* Return value:

void csext_lcg(double defgrad[3][3], double lcg[3][3]);

csext_linsolv: Function to solve a linear system of equations

* Function: csext_linsolv

* -------------------

* Solves a linear system of equations, 'Ax = b', for 'x',

* where 'A' is (n x n), and 'b' is (n x 1).

* The maximum allowed size of the system is n = 6.

* The solution vector 'x' is stored in 'b' as output.

* The elements of matrix 'A' are not preserved.

* If successful, 0 is returned.

* If the solution cannot be determined, -1 is returned.

* If n < 1 or n > 6, -2 is returned.

* double *b (input/output)

* Return value:

CSEXTUTILS_SYMBOLS int csext_linsolv(int n, double *A, double *b);

csext_lpolar: Function that computes the Left polar decomposition of a matrix

* Function: csext_lpolar

* -------------------

* Computes the Left polar decomposition F = VR,

* such that the deformation gradient 'defgrad' is

* multiplicatively decomposed into a rotation, 'R',

* and a stretch tensor, 'V'.

* If the polar decomposition fails, -1 is returned.

* If successful, 0 is returned.

* double defgrad[3][3]

* double V[3][3] (output)

* double R[3][3] (output)

* Return value:

int csext_lpolar(double defgrad[3][3], double V[3][3],

double R[3][3]);

csext_mises: Function that computes the von Mises equivalent stress

* Function: csext_mises

* ---------------------

* Computes and returns the von Mises equivalent stress based

* on a stress tensor 'sig'.

* double sig[3][3]

* Return value:

double csext_mises(double sig[3][3]);

csext_mul: Function that multiplies two matrices

* Function: csext_mul

* -------------------

* Multiplies two (3 x 3) matrices 'A' and 'B'.

* The result 'AB' is stored in 'C'.

* double A[3][3]

* double B[3][3]

* double C[3][3] (output)

* Return value:

void csext_mul(double A[3][3], double B[3][3], double C[3][3]);

csext_rcg: Function that computes the Right Cauchy–Green deformation tensor

* Function: csext_rcg

* -------------------

* Computes the Right Cauchy-Green deformation tensor 'rcg'

* based on the deformation gradient 'defgrad'.

* double defgrad[3][3]

* double rcg[3][3] (output)

* Return value:

void csext_rcg(double defgrad[3][3], double rcg[3][3]);

csext_rpolar: Function that computes the Right polar decomposition of a matrix

* Function: csext_rpolar

* -------------------

* Computes the Right polar decomposition F = RU,

* such that the deformation gradient 'defgrad' is

* multiplicatively decomposed into a rotation, 'R',

* and a stretch tensor, 'U'.

* If the polar decomposition fails, -1 is returned.

* If successful, 0 is returned.

* double defgrad[3][3]

* double R[3][3] (output)

* double U[3][3] (output)

* Return value:

int csext_rpolar(double defgrad[3][3], double R[3][3],

double U[3][3]);

csext_spect: Function that computes a matrix based on its spectral decomposition

* Function: csext_spect

* -------------------

* Computes a symmetric (3 x 3) matrix 'A' based on

* its spectral decomposition A=Q*diag*Q^T.

* The matrix 'diag' is diagonal and stores the

* eigenvalues of 'A'. The matrix 'Q' stores the

* eigenvectors (column-wise) of 'A', with the ordering

* corresponding to the eigenvalues in 'diag'.

* The vector 'd' stores the diagonal elements of 'diag'.

* double Q[3][3]

* double A[3][3] (output)

* Return value:

void csext_spect(double Q[3][3], double d[3], double A[3][3]);

csext_trace: Function that computes the matrix trace

* Function: csext_trace

* ---------------------

* Computes and returns the trace of a (3 x 3) matrix 'A'.

* double A[3][3]

* Return value:

double csext_trace(double A[3][3]);

csext_transp: Function that computes the matrix transpose

* Function: csext_transp

* ----------------------

* Computes the transpose of a (3 x 3) matrix 'A'.

* The result is stored in 'transp'.

* double A[3][3]

* double transp[3][3] (output)

* Return value:

void csext_transp(double A[3][3], double transp[3][3]);

csext_add_voigt: Function that adds two matrices

* Function: csext_add_voigt

* ---------------------------

* Adds two matrices 'A' and 'B' stored on Voigt form.

* The result is stored in 'C'.

* double C[6] (output)

* Return value:

void csext_add_voigt(double A[6], double B[6], double C[6]);

csext_dev_voigt: Function that computes the deviator of a symmetric matrix

* Function: csext_dev_voigt

* ---------------------------

* Computes the deviator of a (symmetric) matrix 'A' stored on Voigt form.

* The result is stored in 'dev'.

* double dev[6] (output)

* Return value:

void csext_dev_voigt(double A[6], double dev[6]);

csext_dot_voigt: Function that computes the inner product of symmetric matrices

* Function: csext_dot_voigt

* ---------------------------

* Computes and returns the dot product (inner product) of two

* symmetric (3 x 3) matrices 'A' and 'B' stored on Voigt form.

* Return value:

double csext_dot_voigt(double A[6], double B[6]);

csext_mises_voigt: Function that computes the von Mises equivalent stress

* Function: csext_mises_voigt

* ---------------------

* Computes and returns the von Mises equivalent stress based

* on a stress tensor 'sig' on Voigt form.

* double sig[6]

* Return value:

double csext_mises_voigt(double sig[6]);

csext_utils_trace_voigt: Function that computes the matrix trace

* Function: csext_trace_voigt

* ---------------------------

* Computes and returns the trace of a symmetric (3 x 3) matrix

* 'A' stored on Voigt form.

* Return value:

double csext_trace_voigt(double A[6]);

csext_from_voigt: Function to change from Voigt notation

* Function: csext_from_voigt

* -------------------

* Converts a symmetric (3 x 3) matrix 'A' stored

* on Voigt form to matrix form. The result is stored in 'B'.

* If 'def' = 1, the values of last three elements of 'A' are

* un-altered when passed into 'B'.

* If 'def' = 2, the values of last three elements of 'A' are

* halved when passed into 'B'.

* Returns -1 if the value of 'def' is invalid, 0 otherwise.

* double B[3][3] (output)

* Return value:

int csext_from_voigt(double A[6], int def, double B[3][3]);

csext_to_voigt: Function to change to Voigt notation

* Function: csext_to_voigt

* -------------------

* Converts a symmetric (3 x 3) matrix 'A' to Voigt form.

* The result is stored in 'B'.

* If 'def' = 1, the values of the off-diagonal components

* of 'A' are un-altered when passed into 'B'.

* If 'def' = 2, the values of the off-diagonal components

* of 'A' are doubled when passed into 'B'.

* Returns -1 if the value of 'def' is invalid, 0 otherwise.

* double A[3][3]

* double B[6] (output)

* Return value:

int csext_to_voigt(double A[3][3], int def, double B[6]);

csext_jac_conv: Function to convert the Jacobian

* Function: csext_jac_conv

* ------------------------

* Converts a Jacobian from 'dSde' to 'dSdF' using the deformation

* gradient 'defgrad', where:

* - 'S' is the 2:nd Piola-Kirchhoff stress tensor,

* - 'e' is the Green-Lagrange strain tensor,

* - 'F' is the deformation gradient.

* If 'def' = 1, dSde is defined using tensor shears.

* If 'def' = 2, dSde is defined using engineering shears.

* Returns -1 if the value of 'def' is invalid.

* double dSde[6][6]

* double defgrad[3][3]

* double dSdF[6][9] (output)

* Return value:

int csext_jac_conv(double dSde[6][6], int def,

double defgrad[3][3], double dSdF[6][9]);

csext_unit: Function to define a unit matrix

* Function: csext_unit

* -------------------

* Initializes an (n x n) matrix 'A' to the identity matrix.

* Returns -1 if n < 1,

* double *A (output)

* Return value:

int csext_unit(int n, double *A);

csext_zero: Function to initialize a matrix

* Function: csext_zero

* -------------------

* Initializes an (m x n) matrix 'A' to zero.

* Returns -1 if m < 1 or n < 1,

* double *A (output)

* Return value:

int csext_zero(int m, int n, double *A);