Constitutive Laws

Every time a ``deformable'' entity requires a constitutive law, a template constitutive law is read. This has been implemented by means of C++ templates in order to allow the definition of a general constitutive law when possible. The ``deformable'' elements at present are:

• rod and genel spring and related elements (1D);
• deformable hinge and deformable displacement joint elements (3D);
• deformable joint and beam elements (6D).

Constitutive laws are also used in non-structural components, to allow some degree of generality in defining input/output relationships. Some constitutive laws are meaningful only when related to some precise dimensionality. In some special cases, general purpose (genel) elements use 1D constitutive laws to express an arbitrary dependence of some value on a scalar state of the system. Table 2.6 shows the availability of each constitutive law.

The meaning of the input and output parameters of a constitutive law is dictated by the entity that uses it. In general, the user should refer to the element the constitutive law is being instantiated for in order to understand what the input and the output parameters are supposed to be.

Usually, constitutive laws can be directly defined when required, according to the definition of an element. However, the special card described in Section 2.3.1 allows to define constitutive laws stand-alone, and attach them to the elements by means of the following mechanism:

    <constitutive_law> ::=
        { <constitutive_law_definition>
            | <constitutive_law_reference> }

    <constitutive_law_reference> ::= reference , <label>

where <constitutive_law_definition> is described in the following, while <label> is the label of a previously defined constitutive law of the appropriate dimensionality, as described in Section 2.3.1.

Table 2.6: Constitutive laws dimensionality

Constitutive law	1D	3D	6D
linear elastic, linear elastic isotropic	$\surd$	$\surd$	$\surd$
linear elastic generic	$\surd$	$\surd$	$\surd$
linear elastic generic axial torsion coupling			$\surd$
cubic elastic generic	$\surd$	$\surd$
log elastic	$\surd$
linear elastic generic bi-stop	$\surd$	$\surd$	$\surd$
double linear elastic	$\surd$	$\surd$
isotropic hardening elastic	$\surd$	$\surd$	$\surd$
scalar function elastic, scalar function elastic isotropic	$\surd$	$\surd$	$\surd$
scalar function elastic orthotropic	$\surd$	$\surd$	$\surd$
linear viscous, linear viscous isotropic	$\surd$	$\surd$	$\surd$
linear viscous generic	$\surd$	$\surd$	$\surd$
linear viscoelastic, linear viscoelastic isotropic	$\surd$	$\surd$	$\surd$
linear viscoelastic generic	$\surd$	$\surd$	$\surd$
linear time variant viscoelastic generic	$\surd$	$\surd$	$\surd$
linear viscoelastic generic axial torsion coupling			$\surd$
cubic viscoelastic generic	$\surd$	$\surd$
double linear viscoelastic	$\surd$	$\surd$
turbulent viscoelastic	$\surd$
linear viscoelastic generic bi-stop	$\surd$	$\surd$	$\surd$
shock absorber	$\surd$
symbolic elastic	$\surd$	$\surd$	$\surd$
symbolic viscous	$\surd$	$\surd$	$\surd$
symbolic viscoelastic	$\surd$	$\surd$	$\surd$
ann elastic	$\surd$	$\surd$	$\surd$
ann viscoelastic	$\surd$	$\surd$	$\surd$
nlsf elastic	$\surd$	$\surd$	$\surd$
nlsf viscous	$\surd$	$\surd$	$\surd$
nlsf viscoelastic	$\surd$	$\surd$	$\surd$
nlp elastic	$\surd$	$\surd$	$\surd$
nlp viscous	$\surd$	$\surd$	$\surd$
nlp viscoelastic	$\surd$	$\surd$	$\surd$
Constitutive law wrapper	1D	3D	6D
invariant angular (wrapper)		$\surd$

The constitutive laws are entered as follows:

    <constitutive_law_definition> ::= <specific_const_law>
        [ , prestress, (entity) <prestress> ]
        [ , prestrain, (entity_tpl_driver) <prestrain> ]
    <specific_const_law> ::= <const_law_name> ,
        <const_law_data>

where const_law_name is the name of the constitutive law and const_law_data depends on the specific constitutive law. The latter fields, whose type depends on the dimensionality of the constitutive law, are optional, under the assumption that the constitutive law is the last portion of a card, or that any ambiguity can be avoided. The data specific to the currently available constitutive laws must be entered as follows:

Linear elastic, linear elastic isotropic

    <specific_const_law> ::= linear elastic [ isotropic ] , 
        (scalar) <stiffness>

the isotropic stiffness coefficient; the word isotropic can be omitted, essentially because it has no meaning for scalar constitutive laws.

Example.

    constitutive law: 1, name, "scalar isotropic law",
        1, linear elastic, 1.e9;
    constitutive law: 2, name, "3D isotropic law",
        3, linear elastic isotropic, 1.e9;
    constitutive law: 3, name, "6D isotropic law",
        6, linear elastic isotropic, 1.e9;

Linear elastic generic

    <specific_const_law> ::= linear elastic generic ,  
        (derivative_of_entity) <stiffness>

the stiffness matrix. In case of 1D, the type is scalar, and there is no distinction between generic and isotropic, while, in case of n x 1 vectors, the type is the corresponding n x n matrix.

Example.

    constitutive law: 1, name, "scalar isotropic law",
        1, linear elastic generic, 1.e9;
    constitutive law: 2, name, "3D isotropic law",
        3, linear elastic generic,
            sym, 1.e9,  0.,    0.,
                        1.e6, -1.e5,
                               1.e6;
    constitutive law: 3, name, "6D isotropic law",
        6, linear elastic generic,
            diag, 1.e9, 1.e9, 1.e9, 1.e6, 1.e6, 1.e6;

Linear elastic generic axial torsion coupling

    <specific_const_law> ::= 
        linear elastic generic axial torsion coupling ,  
            (derivative_of_entity) <stiffness> ,
            (scalar) <coupling_coefficient>

this is defined only for 6 x 1 vectors, where the torsion stiffness, coefficient a₄₄ in the stiffness matrix, depends linearly on the axial strain, $\varepsilon_{1}^{}$, by means of coupling_coefficient, i.e. the run-time torsion stiffness is

a₄₄ = GJ + coupling$$\_coefficient$$ ^. $$\varepsilon_{1}^{}$$.

This <coupling_coefficient>, in the classical nonlinear beam theory, is estimated according to some geometric property [2]; a general approach to the computation of prestressed beam properties is presented in [3], which is implemented in some versions of the ANBA software.

Cubic elastic generic

    <specific_const_law> ::= 
        cubic elastic generic
            (entity) <stiffness1> ,
            (entity) <stiffness2> ,
            (entity) <stiffness3>

this is defined only for scalar and 3 x 1 vectors; the constitutive law is written according to the formula

F = stiffness1 ^. $$\varepsilon$$ + stiffness2 ^. $$\left\vert\vphantom{ \varepsilon }\right.$$$$\varepsilon$$$$\left.\vphantom{ \varepsilon }\right\vert$$$$\varepsilon$$ + stiffness3 ^. $$\varepsilon^{3}_{}$$

and it is mainly intended for use with human body models where the stiffness of the joints is typically given in this form.

Log elastic

    <specific_const_law> ::= log elastic ,
        (derivative_of_entity) <stiffness>

this is defined only for scalars. The force is defined as:

$$f \ = \ \texttt{stiffness} \ \ensuremath{\mathrm{log}} \left( 1+\varepsilon \right)$$

Linear elastic bi-stop generic

    <specific_const_law> ::= linear elastic bistop,
        (derivative_of_entity) <stiffness> ,
        [ initial state , { inactive | active } , ]
        (DriveCaller)<activating_condition> ,
        (DriveCaller)<deactivating_condition>

Double linear elastic

    <specific_const_law> ::= double linear elastic ,
        (scalar) <stiffness_1> ,
        (scalar) <upper_strain> ,
        (scalar) <lower_strain> ,
        (scalar) <stiffness_2>

this is defined for scalar and 3 x 1 vectors. In the scalar case the meaning of the entries is straightforward, while in case of 3 x 1 vectors, the constitutive law is isotropic but in the local direction 3, where, in case of strain out of the upper or lower bound, the stiffness_2 is used.

Isotropic hardening elastic

    <specific_const_law> ::= isotropic hardening elastic ,
        (scalar) <stiffness> ,
        (scalar) <reference_strain>
        [ , linear stiffness , <linear_stiffness> ]

this constitutive law is defined as follows:

f  =  stiffness $${\frac{{ \beta + \alpha \left\vert \boldsymbol{\varepsilon} \rig... ...rt ^2 }}{{ 1 + \alpha \left\vert \boldsymbol{\varepsilon} \right\vert ^2 }}}$$$$\varepsilon$$

where $\alpha$ = 3/$\left\vert\vphantom{ \texttt{reference\_strain} }\right.$reference$\_strain$$\left.\vphantom{ \texttt{reference\_strain} }\right\vert^{2}_{}$, and $\beta$ = linear$\_stiffness$/stiffness. The resulting constitutive law, in the scalar case, is somewhat soft/hard when $\beta$ is lower/greater than 1 and $\left\vert\vphantom{ \varepsilon }\right.$$\varepsilon$$\left.\vphantom{ \varepsilon }\right\vert$ is smaller than reference_strain, approaching linear_stiffness when $\varepsilon$ $\rightarrow$ 0, while it grows to quasi-linear for larger $\left\vert\vphantom{ \varepsilon }\right.$$\varepsilon$$\left.\vphantom{ \varepsilon }\right\vert$, with slope stiffness.

Scalar function elastic, scalar function elastic isotropic

This constitutive law is based on a scalar function (Section ) to represent an analytical force-displacement curve of a single variable that is automatically differentiated to compute the slope of the curve, namely the local stiffness.

    <specific_const_law> ::= scalar function elastic isotropic ,
        (DifferentiableScalarFunction)<function>

this constitutive law is defined as follows:

f_i = < function >$$\left(\vphantom{ \boldsymbol{\varepsilon}_i }\right.$$$$\boldsymbol\varepsilon$$_i$$\left.\vphantom{ \boldsymbol{\varepsilon}_i }\right)$$

the force is computed for each direction as a function of the respective strain component using the same function. When used for 1D elements, the word isotropic can be omitted.

Scalar function elastic orthotropic

    <specific_const_law> ::= scalar function elastic orthotropic ,
        { (DifferentiableScalarFunction)<function> | null }
        [ , ... ]

this constitutive law is defined as follows:

f_i = < function >_i$$\left(\vphantom{ \boldsymbol{\varepsilon}_i }\right.$$$$\boldsymbol\varepsilon$$_i$$\left.\vphantom{ \boldsymbol{\varepsilon}_i }\right)$$

the force is computed for each direction as a function of the respective strain component using a specific function for each component. If no force ought to be used for a direction, for example because that direction is constrained by a kinematic joint, the keyword null can be used to indicate that no function is expected for that component. When used for 1D elements, the word orthotropic can be omitted; note that in this case a scalar function elastic isotropic constitutive law is actually instantiated.

Example.

    scalar function: "myfunc", multilinear,
        -1., -100.,
        -.5, -70.,
        0., 0.,
        .5, 70.,
        1., 100.;
    constitutive law: 1000, 3,
        scalar function elastic orthotropic,
            null,
            null,
            "myfunc";

indicates that the constitutive law is only defined in direction 3 as a multilinear function.

Linear viscous, linear viscous isotropic

    <specific_const_law> ::= linear viscous [ isotropic ] , 
        (scalar) <viscosity_coefficient>

the linear viscous coefficient.
Note: this constitutive law does not require any prestrain template drive caller.

Linear viscous generic

    <specific_const_law> ::= linear viscous generic , 
        (derivative_of_entity) <viscosity_matrix>

the linear viscous matrix.
Note: this constitutive law does not require any prestrain template drive caller.

Linear viscoelastic, linear viscoelastic isotropic

    <specific_const_law> ::= linear viscoelastic [ isotropic ] ,
        (scalar) <stiffness> ,
        { (scalar) <viscosity_coefficient>
            | proportional, (scalar) <factor> }

the isotropic stiffness and viscosity coefficients.

Linear viscoelastic generic

    <specific_const_law> ::= linear viscoelastic generic ,  
        (derivative_of_entity) <stiffness> ,
        { (derivative_of_entity) <viscosity_matrix> 
            | proportional, (scalar) <factor> }

the linear stiffness and viscosity matrices.

Linear time variant viscoelastic generic

    <specific_const_law> ::= linear time variant viscoelastic generic ,  
        (derivative_of_entity) <stiffness> ,
        (DriveCaller> <stiffness_scale> ,
        { (derivative_of_entity) <viscosity>
            | proportional, (scalar) <factor> } ,
        (DriveCaller> <viscosity_scale>

the linear stiffness and viscosity matrices are multiplied by the respective scale factors,

$$\_scale \boldsymbol\varepsilon \_scale \dot{{\boldsymbol{\varepsilon}}}$$ f (2.40)

Example.

    linear time variant visco elastic,
        1000., cosine, 2., pi/.2, .1/2, half, 1.,
        100., cosine, 2., pi/.2, 1/2, half, 1.

At 2s, the stiffness grows of 10% from the nominal value in .2s, while the damping doubles. See Section 2.5 for details on the syntax of drive callers.

Beware that arbitrarily changing the stiffness and the damping of an elastic component during the execution of the simulation may have no physical meaning. The intended use of this feature is for tailoring the analysis; for example, a higher damping level may be desirable to smooth out a non-physical transient, and later return to the appropriate damping value.

Linear viscoelastic generic axial torsion coupling

    <specific_const_law> ::= 
        linear viscoelastic generic axial torsion coupling ,  
            (derivative_of_entity) <stiffness> ,
            { (derivative_of_entity) <viscosity_matrix> 
                | proportional, (scalar) <factor> }
            (scalar) <coupling_coefficient>

this is defined only for 6 x 1 vectors; it is the viscoelastic extension of the linear elastic generic axial torsion coupling constitutive law.

Cubic viscoelastic generic

    <specific_const_law> ::= 
        cubic elastic generic
            (entity) <stiffness1> ,
            (entity) <stiffness2> ,
            (entity) <stiffness3> ,
            (derivative_of_entity) <viscosity_matrix>

this is defined only for scalar and 3 x 1 vectors; the constitutive law is written according to the formula

F = stiffness1 ^. $$\varepsilon$$ + stiffness2 ^. $$\left\vert\vphantom{ \varepsilon }\right.$$$$\varepsilon$$$$\left.\vphantom{ \varepsilon }\right\vert$$$$\varepsilon$$ + stiffness3 ^. $$\varepsilon^{3}_{}$$ + viscosity$$\_matrix$$ ^. $$\dot{{\varepsilon}}$$

and it is mainly intended for use with human body models where the stiffness of the joints is typically given in this form.

Double linear viscoelastic

    <specific_const_law> ::= double linear viscoelastic ,
        (scalar) <stiffness_1> ,
        (scalar) <upper_strain> ,
        (scalar) <lower_strain> ,
        (scalar) <stiffness_2> ,
        (scalar) <viscosity_coefficient>
        [ , second damping , (scalar) <viscosity_coefficient_2> ]

this is analogous to the double linear elastic constitutive law, except for the isotropic viscosity term. The second viscosity value is used when the strain is outside the lower_strain-upper_strain range.

When this constitutive law is used with 3 x 1 vectors, the double linear elastic and viscous term only applies to component 3.

Turbulent viscoelastic

    <specific_const_law> ::= turbulent viscoelastic ,
        (scalar) <stiffness> ,
        (scalar) <parabolic_viscous_coefficient>
        [ , (scalar) <threshold> 
            [ , (scalar) <linear_viscous_coefficient> ] ]

the constitutive law has the form:

f  =  stiffness $$\boldsymbol\varepsilon$$ + k $$\dot{{\boldsymbol{\varepsilon}}}$$

where:

k = $$\left\{\vphantom{ \begin{array}{ lcr } \texttt{linear\_viscous\... ...ymbol{\dot{\varepsilon}} \right\vert > \texttt{threshold} \end{array} }\right.$$$$\begin{array}{ lcr } \texttt{linear\_viscous\_coefficient} & & ... ...rt \boldsymbol{\dot{\varepsilon}} \right\vert > \texttt{threshold} \end{array}$$

if threshold is null, or not defined, the constitutive law is always parabolic. If the linear_viscous_coefficient is not defined, it is computed based on parabolic_viscous_coefficient and on threshold to give a continuous force curve (with discontinuous slope). Otherwise, it can be set by the user to give a discontinuous force curve, as observed in some fluids at intermediate Reynolds number.

Linear viscoelastic bi-stop generic

    <specific_const_law> ::= linear viscoelastic bistop ,
        (derivative_of_entity) <stiffness> ,
        (derivative_of_entity) <viscosity_coefficient> ,
        [ initial state , { inactive | active } , ]
        (DriveCaller)<activating_condition> ,
        (DriveCaller)<deactivating_condition>

GRAALL damper

This is a very experimental constitutive law, based on a nonlinear model for a hydraulic damper to be used in landing gear modeling. Basically, it requires the user to supply the name of the GRAALL-style input file with damper data. It will be documented as soon as it reaches an appreciable level of stability. See also the shock absorber constitutive law.

shock absorber

This constitutive law implements a landing gear hydraulic shock absorber:

    <specific_const_law> ::= shock absorber ,
        [ prestrain , <value> , ]
        <reference pressure> ,
        <reference area for force computation> ,
        <interaction coefficient> ,
        <polytropic exponent> ,
        [ epsilon max , <upper strain bound> , ]
        [ epsilon min , <lower strain bound> , ]
        [ penalty , <penalty factor for strain> , 
            <penalty factor for strain rate> , ]
        [ metering , <metering area> ,
            [ negative , <metering area for negative strain rate> , ]
        [ orifice , <orifice area> , ]
        <fluid area> ,
        <fluid density> ,
        <drag coefficient / reference length
                (scales strain rate to velocity)>
        [ , friction, <reference epsilon prime> ,
                 <friction amplitude coefficient> ]

where

• the interaction coefficient is represented by
  kinematicscale$${\frac{{L A}}{{V_0}}}$$
  where kinematic scale is the ratio between the stroke of the shock absorber and that of the gas;
• epsilon max is the upper strain bound; it must be at least larger than the prestrain, and defaults to 0, i.e. the shock absorber, at rest, is assumed to be fully extended;
• epsilon min is the lower strain bound; it must be at least smaller than the prestrain, and defaults to -0.5, i.e. the shock absorber is assumed to allow a contraction equal to half its full length;
• the penalty factor for strain defaults to 1e+9; it is active only when strain bounds are violated;
• the penalty factor for strain rate defaults to 0; it is active only when strain bounds are violated;
• the metering area is given by a drive and is strain dependent; if the keyword negative is used, then the metering area for negative strain rate is used when the strain rate is negative, i.e. the shock absorber is being compressed, while the metering area is used only when the shock absorber is extending;
• the <orifice> drive determines the area of an additional orifice, which essentially depends on the sign of the strain rate; it is used to implement relief valves;
• ...

This constitutive law contributes to the output of the element it is associated with. The name associated to each output contribution can be used to reference the corresponding value as element private data. The contributions are

1. "p" gas pressure;
2. "A" metering area;
3. "Fe" elastic force;
4. "Fv" viscous force.

symbolic elastic

The implementation of the family of symbolic constitutive laws is based on GiNaC (http://www.ginac.de/), a free software package for symbolic algebra manipulation. It is essentially used to automatically differentiate the user-supplied expression that describes the relationship between the output and the input.

The syntax is

    <specific_const_law> ::= symbolic elastic ,
        epsilon , " <epsilon> " ,
        expression , " <expression> "

where epsilon is the symbol describing the input parameter as it will be used in expression.

For constitutive laws with more than one dimension, a string for each epsilon and one for each expression are expected. For example:

    # 1D symbolic constitutive law
    constitutive law: 1001, 1, symbolic elastic,
        epsilon, "eps",
        expression, "1000.*eps + 5.*eps^3";
    # 3D symbolic constitutive law
    constitutive law: 1003, 3, symbolic elastic,
        epsilon, "eps1", "eps2", "eps3",
        expression,
            "1000.*eps1 + 5.*eps1^3 - 10.*eps2*eps3",
            "1000.*eps2 + 5.*eps2^3 - 10.*eps3*eps1",
            "1000.*eps3 + 5.*eps3^3 - 10.*eps1*eps2";

Note: right now, the symbols defined within the mathematical parser are not available within symbolic consitutive laws.

symbolic viscous

The syntax is

    <specific_const_law> ::= symbolic viscous ,
        epsilon prime , " <epsilon_prime> " ,
        expression , " <expression> "

where epsilon_prime is the symbol describing the derivative of the input parameter as it will be used in expression.

For constitutive laws with more than one dimension, a string for each epsilon_prime and one for each expression are expected.

symbolic viscoelastic

The syntax is

    <specific_const_law> ::= symbolic viscoelastic,
        epsilon , " <epsilon> " ,
        epsilon prime , " <epsilon_prime> " ,
        expression , " <expression> "

where epsilon and epsilon_prime are the symbols describing the input parameter and its derivative as they will be used in expression.

For constitutive laws with more than one dimension, a string for each epsilon, one for each epsilon_prime, and one for each expression are expected.

ann elastic

The implementation of the family of ann constitutive laws is based on an Artifical Neural Network library that is embedded in MBDyn.

The syntax is

    <specific_const_law> ::= ann elastic ,
        " <file_name> "

where the file file_name contains the parameters of the network. What is mainly significant to users is the need to scale inputs and outputs to match the amplitudes and possibly the offsets used in the training. For this purpose, the last two rows of the input files contain coefficients

    b1 b0
    a1 a0

which are used to scale the input u and the output y according to the transformation

$$\overline{{u}}$$ = b₁u + b₀ (2.41)

$$\overline{{y}}$$ = a₁y + a₀ (2.42)

so that the actual output is

$${\frac{{f \left( b_1 u + b_0 \right) - a_0}}{{a_1}}}$$ (2.43)

ann viscoelastic

The syntax is

    <specific_const_law> ::= ann viscoelastic ,
        " <file_name> "

where the file file_name contains the parameters of the network.

nlsf viscoelastic

This constitutive law was sponsored by Hutchinson CdR.

For an n-dimensional constitutive law, it implements the formula

$$\sum_{{j=1,n}}^{} \varepsilon_{j}^{} \left(\vphantom{ \varepsilon_i }\right. \varepsilon_{i}^{} \left.\vphantom{ \varepsilon_i }\right) \sum_{{j=1,n}}^{} \dot{{\varepsilon}}_{j}^{} \left(\vphantom{ \dot{\varepsilon_i} }\right. \dot{{\varepsilon_i}} \left.\vphantom{ \dot{\varepsilon_i} }\right)$$ (2.44)

where f'_i and f''_i are arbitrary instances of the scalar functions primitive, which includes piecewise linear and spline regularization of given data, while k'_0ij and k''_0ij are the constant coefficients of a linear viscoelastic model, that account for the cross-couplings between the stresses and the strains and strain rates.

Elastic and viscous variants are defined. They differ from the viscoelastic one by only allowing the specific fraction of the input data.

Syntax.

The syntax is

    <specific_const_law> ::= nlsf viscoelastic ,
        (derivative_of_entity) <kappa_0'> ,
        { null | (scalar_function) <diag_force'> }
        [ , ... ] ,
        { (derivative_of_entity) <kappa_0''> | proportional , <coef> } ,
        { null | (scalar_function) <diag_force''> }
        [ , ... ]

    <specific_const_law> ::= nlsf elastic ,
        (derivative_of_entity) <kappa_0'> ,
        { null | (scalar_function) <diag_force'> }
        [ , ... ]

    <specific_const_law> ::= nlsf viscous ,
        (derivative_of_entity) <kappa_0''> ,
        { null | (scalar_function) <diag_force''> }
        [ , ... ]

The terms kappa_0' and kappa_0'' are indicated as ``derivative of entity'' because in the generic case they are the result of a differential operator that computes the derivative of the force vector as function of the strain or strain rate vectors. So, for a 3D constitutive law, the force, the strain and the strain rate are 3 x 1 vectors, while kappa_0' and kappa_0'' are 3 x 3 matrices. Matrix input, in MBDyn, requires to write the whole set of coefficients, row-wise; however, it can be described in many synthetic manners if matrices have some special properties, as typical linear constitutive laws do.

For example, symmetric, diagonal and empty matrices have specific short forms; in the 3 x 3 case they are

        # generic 3x3 matrix:
        10., -2., -2.,
        -2., 20., -8.,
        -2., -8., 20.

        # symmetric 3x3 matrix (same as above matrix):
        sym,
        10., -2., -2.,
             20., -8.,
                  20.

        # diagonal 3x3 matrix (diagonal of above matrix):
        diag,
        10., 20., 20.

        # empty 3x3 matrix:
        null

The terms diag_force' and diag_force'' refer to the names of the scalar functions that are defined for each component of the force vector. There must be as many terms as the dimensions of the constitutive law. These terms can either be null, if no non-linear contribution is defined for that component, or contain a string, enclosed in double quotes, which must refer to an already defined scalar function. The same scalar function can be used multiple times.

Example.

A 1D and a 3D constitutive law using the same exponential scalar function:

    # define a scalar function
    scalar function: "exponential", exp, coefficient, -2., 5.;

    # define a 1D constitutive law
    constitutive law: 1000, 1,
        nlsf viscoelastic,
            10.,
                null,               # stiffness is just linear
            0.,
                "exponential";      # damping is exponential

    # define a 3D constitutive law
    constitutive law: 3000, 3,
        nlsf viscoelastic,
            sym,
                10., -2., -2.,
                     10., -2.,
                          10.,
                null,               # stiffness is just linear
                null,               # stiffness is just linear
                null,               # stiffness is just linear
            null,
                "exponential",      # damping is exponential...
                null,               # ...but in direction 1 only!
                null;

The first constitutive law corresponds to

$$\varepsilon \dot{{\varepsilon}}$$ (2.45)

while the second one corresponds to

$$\left[\vphantom{ \begin{array}{ ccc } 10.0 & -2.0 & -2.0 \\ -2.0 & 10.0 & -2.0 \\ -2.0 & -2.0 & 10.0 \end{array} }\right. \begin{array}{ ccc } 10.0 & -2.0 & -2.0 \\ -2.0 & 10.0 & -2.0 \\ -2.0 & -2.0 & 10.0 \end{array} \left.\vphantom{ \begin{array}{ ccc } 10.0 & -2.0 & -2.0 \\ -2.0 & 10.0 & -2.0 \\ -2.0 & -2.0 & 10.0 \end{array} }\right] \left\{\vphantom{ \begin{array}{c} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3 \end{array} }\right. \begin{array}{c} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3 \end{array} \left.\vphantom{ \begin{array}{c} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3 \end{array} }\right\} \left\{\vphantom{ \begin{array}{c} 5.0 \cdot \mathrm{e}^{-2 \dot{\varepsilon}_1} \\ 0.0 \\ 0.0 \end{array} }\right. \begin{array}{c} 5.0 \cdot \mathrm{e}^{-2 \dot{\varepsilon}_1} \\ 0.0 \\ 0.0 \end{array} \left.\vphantom{ \begin{array}{c} 5.0 \cdot \mathrm{e}^{-2 \dot{\varepsilon}_1} \\ 0.0 \\ 0.0 \end{array} }\right\}$$ (2.46)

nlp viscoelastic

This constitutive law was sponsored by Hutchinson CdR.

This constitutive law consists in the linear combination of nonlinear elastic and viscous effects whose coefficients are arbitrary nonlinear functions of the strain $\varepsilon$:

$$\sum_{{j=1,3}}^{} \left(\vphantom{ k'_{0ij} + \delta_{ij} k_i' \left( \varepsilon_i \right) }\right. \delta_{{ij}}^{} \left(\vphantom{ \varepsilon_i }\right. \varepsilon_{i}^{} \left.\vphantom{ \varepsilon_i }\right) \left.\vphantom{ k'_{0ij} + \delta_{ij} k_i' \left( \varepsilon_i \right) }\right) \varepsilon_{j}^{} \sum_{{j=1,3}}^{} \left(\vphantom{ k''_{0ij} + \delta_{ij} k''_i \left( \varepsilon_i \right) }\right. \delta_{{ij}}^{} \left(\vphantom{ \varepsilon_i }\right. \varepsilon_{i}^{} \left.\vphantom{ \varepsilon_i }\right) \left.\vphantom{ k''_{0ij} + \delta_{ij} k''_i \left( \varepsilon_i \right) }\right) \dot{{\varepsilon}}_{j}^{}$$ (2.47)

where $\delta_{{ij}}^{}$ is Kronecker's operator; k'_i and k''_i are implemented as arbitrary scalar functions, which includes piecewise linear and spline regularization of given data.

Elastic and viscous variants are defined. They differ from the viscoelastic one by only allowing the specific fraction of the input data.

Syntax.

The syntax is

    <specific_const_law> ::= nlp viscoelastic ,
        (derivative_of_entity) <kappa_0'> ,
        { null | (scalar_function) <diag_stiffness> }
        [ , ... ] ,
        { (derivative_of_entity) <kappa_0''> | proportional , <coef> } ,
        { null | (scalar_function) <diag_damping> }
        [ , ... ]

    <specific_const_law> ::= nlp elastic ,
        (derivative_of_entity) <kappa_0'> ,
        { null | (scalar_function) <diag_stiffness> }
        [ , ... ]

    <specific_const_law> ::= nlp viscous ,
        (derivative_of_entity) <kappa_0''> ,
        { null | (scalar_function) <diag_damping> }
        [ , ... ]

Example.

A 1D and a 3D constitutive law using the same exponential scalar function:

    # define a scalar function
    scalar function: "exponential", exp, coefficient, -2., 5.;

    # define a 1D constitutive law
    constitutive law: 1000, 1,
        nlp viscoelastic,
            10.,
                null,               # stiffness is just linear
            0.,
                "exponential";      # damping slope is exp.

    # define a 3D constitutive law
    constitutive law: 3000, 3,
        nlp viscoelastic,
            sym,
                10., -2., -2.,
                     10., -2.,
                          10.,
                null,               # stiffness is just linear
                null,               # stiffness is just linear
                null,               # stiffness is just linear
            null,
                "exponential",      # damping slope is exp...
                null,               # ...but in direction 1 only!
                null;

The first constitutive law corresponds to

$$\varepsilon \varepsilon \dot{{\varepsilon}}$$ (2.48)

while the second one corresponds to

$$\left[\vphantom{ \begin{array}{ ccc } 10.0 & -2.0 & -2.0 \\ -2.0 & 10.0 & -2.0 \\ -2.0 & -2.0 & 10.0 \end{array} }\right. \begin{array}{ ccc } 10.0 & -2.0 & -2.0 \\ -2.0 & 10.0 & -2.0 \\ -2.0 & -2.0 & 10.0 \end{array} \left.\vphantom{ \begin{array}{ ccc } 10.0 & -2.0 & -2.0 \\ -2.0 & 10.0 & -2.0 \\ -2.0 & -2.0 & 10.0 \end{array} }\right] \left\{\vphantom{ \begin{array}{c} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3 \end{array} }\right. \begin{array}{c} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3 \end{array} \left.\vphantom{ \begin{array}{c} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3 \end{array} }\right\} \left[\vphantom{ \begin{array}{ ccc } 5.0 \cdot \mathrm{e}^{-2 \v... ...lon_1} & 0.0 & 0.0 \\ 0.0 & 0.0 & 0.0 \\ 0.0 & 0.0 & 0.0 \end{array} }\right. \begin{array}{ ccc } 5.0 \cdot \mathrm{e}^{-2 \varepsilon_1} & 0.0 & 0.0 \\ 0.0 & 0.0 & 0.0 \\ 0.0 & 0.0 & 0.0 \end{array} \left.\vphantom{ \begin{array}{ ccc } 5.0 \cdot \mathrm{e}^{-2 \v... ...lon_1} & 0.0 & 0.0 \\ 0.0 & 0.0 & 0.0 \\ 0.0 & 0.0 & 0.0 \end{array} }\right] \left\{\vphantom{ \begin{array}{c} \dot{\varepsilon}_1 \\ \dot{\varepsilon}_2 \\ \dot{\varepsilon}_3 \end{array} }\right. \begin{array}{c} \dot{\varepsilon}_1 \\ \dot{\varepsilon}_2 \\ \dot{\varepsilon}_3 \end{array} \left.\vphantom{ \begin{array}{c} \dot{\varepsilon}_1 \\ \dot{\varepsilon}_2 \\ \dot{\varepsilon}_3 \end{array} }\right\}$$ (2.49)

Note: although the syntax of this constitutive law and that of the nlsf viscoelastic one (Section 2.9.28) is essentially identical, their behavior may be quite different, as indicated by Equations (2.44) and (2.47).

invariant angular

This is not a constitutive law, but rather a wrapper for constitutive laws used within the ``attached'' variant of the deformable hinge joint. As such, it can only be used with the 3D dimensionality. It basically allows to refer an ancillary 3D constitutive law to an orientation that is intermediate with respect to the orientations of the two connected nodes.

    <specific_const_law> ::= invariant angular ,
        <xi> , <ancillary_const_law>

xi is the fraction of relative orientation the ancillary constitutive law will be referred to, with respect to the first node of the joint. Its value should be comprised between 0 (attached to the first node) and 1 (attached to the second node); a value of 1/2 yields the invariant deformable hinge joint, but any value is allowed, not limited to within the [0, 1] range.

Note: the contribution to the Jacobian matrix will be incomplete when the underlying constitutive law has a viscous contribution. This may slow down convergence, or even prevent it. If this is an issue, the invariant deformable hinge should be used, since it does not suffer from this limitation, although it only allows xi = 1/2.