Materials

The external functions may be written in Fortran, C or C++. The resulting libraries and files should be accessible by OptiStruct regardless of the coding language, provided that consistent function prototyping is respected, and adequate compiling and linking options are used.

Write External Functions

Two Fortran subroutines are required to define user material in OptiStruct. First, a Nonlinear subroutine for the nonlinear large displacement solution, and another subroutine for the linear and nonlinear small displacement solution. Both subroutines are mandatory, and the same argument order should be followed as:

      subroutine usermaterial(idu,stress,
     $     strain,dstrain,dfgrOld,dfgrNew,
     $     stater,state,nstate,
     $     drot,
     $     props,nprops,ndi,nshear,ntens,temp,dtemp,
     $     ieuid,
     $     kinc,dt,t_step,t_total,
     $     cdev,cbulk)

integer idu, nstate, nprops, kinc, ieuid, ndi, nshear, ntens, ierr
double precision stress(6),stater(*),state(*),
     $     cdev(6,6),cbulk, drot(3,3), temp, dtemp, dt, 
     $     t_step, t_total,
     $     strain(6), dstrain(6), props(nprops)
     $     dfgrOld(3,3),dfgrNew(3,3)
character*32000 userdata

subroutine smatusr(idu, nprop, prop, ndi, nshear, ntens, smat, userdata,  ierr)
integer idu, nprop, ndi, nshear, ntens, ierr
double precision prop(nprop), smat(*)
character*32000 userdata

subroutine initusr(idu,nstate,cstate)
      implicit none
      integer idu,nstate
      character*64 cstate(nstate)

      character*64 cstate(nstate)

	cstate(1) = “Variable Name for State(1)”

Subroutine Arguments

Build External Libraries for User-defined Materials

Allows building shared libraries on Windows or Linux.

Refer to Build External Libraries for more information.

Write External Functions

subroutine usrmatht(idu, matID, ieuid, 
				Integration_Point, Layer, Simpson_Point,	
Stater, State,
 nState, drot, props, nprops, temp, 
 dtemp, kinc, dt, t_step, t_total, hgen, smat
 userdata, ierr)

      character*32000 userdata

integer*8 i, idu, matID, ieuid, nState, nprops, kinc, ierr
integer Integration_Point, Layer, Simpson_Point
      double precision Stater(nState), State(nState)
      double precision drot(3,3), rmat(nprops)
      double precision etempmat, dtemp, dt, t_step, t_total
      double precision hgen, smat(9)

subroutine initusr(idu,nstate,cstate)
      implicit none
      integer idu,nstate
      character*64 cstate(nstate)

      character*64 cstate(nstate)

	cstate(1) = “Variable Name for State(1)”

Subroutine Arguments

smat(1): KXX Thermal Conductivity
smat(2): KXY Thermal Conductivity
smat(3): KXZ Thermal Conductivity
smat(4): KYY Thermal Conductivity
smat(5): KYZ Thermal Conductivity
smat(6): KZZ Thermal Conductivity
smat(7): CP Heat Capacity per unit mass
smat(8): RHO Density
smat(9): Free Convection Heat Transfer 
Coefficient (Unused).

Output

Regular Built-in OptiStruct output, like Grid Temperatures, Element Fluxes, and so on can be requested for any thermal material user-subroutine based model. Additionally, user-defined output can also be requested via the State(*) variable in the user-subroutine. The number of such State(*) variables should be identified on the NDEPVAR field.

This State variable can be assigned to any output variable that is required to be output from the subroutine, and it can then be visualized in HyperView after loading the H3D file.

Build External Libraries for User-defined Materials

You can build shared libraries on Windows or Linux.

Refer to Build External Libraries for more information.

It consists of two nonlinear hardening rules, the nonlinear kinematic (NLKIN) and nonlinear isotropic (NLISO) hardening methods.

Generally, the isotropic part is closely related to the von Mises criteria, and the kinematic part is described by the evolution law of back stress.

Combined hardening can be activated by setting HR=6 on the MATS1 Bulk Data.

Isotropic Hardening (NLISO): Nonlinear Yield Function

The yield function of von Mises plasticity can be expressed in a general form as:

The flow rule is defined as change of plastic strain, expressed in rate form as:

Where, $\lambda$ is the rate of plastic multiplier, which is also the rate of equivalent plastic strain.

The flow direction, $N$ can be introduced which is the derivative of the yield function with respect to the stress tensor,

Where, $\eta$ is the relative stress tensor, which is the difference between deviatoric stress and back stress.

For the nonlinear isotropic hardening, the yield stress is assumed to be a power law function of equivalent plastic strain:

Where, $Q$ and $b$ are two parameters which can be directly input via the Q and B fields on the MATS1 data (HR=6, TYPISO=PARAM) or these parameters are computed by parameter fitting algorithms, based on the stress-strain curve from experiment. The isotropic part of the yield stress and the equivalent plastic stress are provided via the SIG and EPS fields on the MATS1 data (HR=6, TYPISO=TABLE).

Nonlinear isotropic hardening is based on the von Mises plasticity criteria and; therefore, is associated with the flow rule.

Kinematic Hardening (NLKIN): Evolution Law of Back Stress

Compared with traditional linear hardening (HR=1 or 2), or the mixed hardening (HR=3), the main difference of NLKIN via HR=6 is the extended evolution law of back stress, which consists of a set of evolution equations for each back stress component:

and

As the evolution of back stress ${\dot{\alpha }}_k$ depends both on the flow direction that is parallel to ${\dot{\epsilon }}_p$ and the back stress component ${\alpha }_k$ itself; thus, the evolution law of back stress is non-associated. This leads to the unsymmetric elasto-plastic consistent tangent modulus. The unsymmetric solver is always turned on as long as combined hardening is active.

Time Integration Scheme

The plasticity problem is usually solved using the return mapping method, which means it is first assumed to be an elastic trial stage, and then the stress is pulled back onto the yield surface if plastic flow occurs. Backward-Euler algorithm is used during the return mapping process.

Parameter Fitting

If the stress-strain curve is provided (TYPKIN=HALFCYCL or TYPISO=TABLE), the data points provided are used to compute the optimal parameters for combined hardening. For parameter fitting, the Levenberg-Marquardt method is used, which is an extension of Newton method. The fitted parameters are printed in the .out file.

For NLISO with TYPISO=TABLE, the provided data in the continuation line is the yield stress versus equivalent plastic strain. This curve is usually generated based on the cyclic experiment with constant strain. The same curve is used for isotropic hardening (HR=1) for example, with TABLES1/TABLEG or TABLEST input.

For NLKIN with TYPKIN=HALFCYCL, the provided data is the yield stress versus equivalent plastic strain, where the yield stress is the total yield stress that is measured directly from experiment, and the equivalent plastic strain is modified by subtracting the elastic strain. For computing the parameters Ci and Gi, the isotropic part will be first subtracted from the curve, and the difference is the hardening part due to kinematic hardening. The subtracted data will be used for parameter fitting. Thus, the provided yield stress values for NLKIN with TYPKIN=HALFCYCL, should always larger than those for NLISO with TYPISO=TABLE.

Temperature-dependent combined hardening

If temperature-dependent combined hardening is active, then all the parameters are temperature dependent, for example,

$\begin{array}{l}{C}_k(T)\\ {\gamma }_k(T)\\ {\sigma }_{y0}(T)\\ Q(T)\\ b(T)\end{array}$

From experiment is only possible to test a limited number of temperatures. Interpolation is used to solving for plasticity at a temperature in between the test temperatures.

Using isotropic hardening NLISO as an example, if the parameters are provided at two temperatures $T_1$ and $T_2$ , then there are two yield stress functions for interpolation:

The interpolated yield stress at the current temperature $T_{curr}$ in [ $T_1$ , $T_2$ ] is then,

with two factors for different temperatures $f_1=\frac{T_2-T_{curr}}{T_2-T_1}$ , $f_2=\frac{T_{curr}-T_1}{T_2-T_1}$ .

As an example, Figure 2 illustrates the principle of interpolation.

It is worthy to mention that the parameters (Q and B) are not interpolated directly. Instead, the yield functions at two temperatures are interpolated. This is the similar case for NLKIN, where the evolution laws at two temperatures are interpolated.

If the temperature, $T$ is beyond the available temperature range [ $T_1$ , $T_2$ ], then closest temperature will be selected. In other words, no extrapolation of temperature is performed. Furthermore, it is suggested to provide the material parameter of NLKIN and NLISO at the same temperature.

In this material model, rate form is used for strain evolution and the strain rate is decomposed into two parts (elastic and plastic parts) in an additive way.

Yield Function

The yielding criteria is defined as:

Where,

$f=\{\begin{array}{c}\frac{2}{3}q\cos \theta +p-{\sigma }_{t}p>-\frac{{\sigma }_c}{3}\\ q-{\sigma }_{c}p\le -\frac{{\sigma }_c}{3}\end{array}$

The first criterion is Rankine criterion (max. Principal stress criterion) for tensile yielding, and the second criterion is von Mises criterion for compressive yielding. Together they form a composite yield surface.

Since the material is assumed to be isotropic, the yield surface can be expressed with a function of three invariant measure of the stress tensor, $p$ , $q$ , and $r$ .

Where, $q$ is von Mises stress of a given stress state.

Where, $p$ is mean stress given by:

The invariant $r$ can be expressed as:

Where, $S$ is the tensor form of deviatoric stress.

The quantity $\theta$ is Lode’s Angle that is defined by $q$ and $r$ as shown below, which identifies the meridional plane for a given stress state.

In the principal stress space, the yield function is a von Mises cylinder intersected with Rankine planes (Figure 3). The yield surface is illustrated in $\pi$ plane in Figure 4 and Figure 5, and in meridional plane in Figure 6.

Where, $S$ is the tensor form of deviatoric stress.

Yield Surface of cast iron plasticity.

Flow Rule

The flow rule is expressed as:

Hardening Rule

In this material model, only isotropic hardening is considered. The yield strength is increasing as the equivalent plastic strain grows. It is noted that as tension and the compression have different flow rules, the equivalent plastic strain is different for tension and compression. Hence, the tensile yield strength and the compressive yield strength are looked up from the tensile and the compressive table separately.

The compressive equivalent plastic strain is expressed as:

Where, $S\left({\dot{\epsilon }}^p\right)$ is the tensor of deviatoric part of plastic strain rate.

The tensile equivalent plastic strain is:

Where, ${\dot{\epsilon }}_{vol}^p={\dot{\epsilon }}_{xx}^p+{\dot{\epsilon }}_{yy}^p+{\dot{\epsilon }}_{zz}^p$ .

In the above expressions, the subscript index $n$ implies the physical quantity is at time $t=n$ , and $n-1$ means the physical quantity is at the last converged time $t=n-1$ .

Input

Cast iron plasticity material data can be specified by using the MCIRON Bulk Data Entry which has the same identification number as an existing MAT1 Bulk Data Entry. The Tensile and Compressive stress-strain curves can be defined using the TABLE_T and TABLE_C fields on the MCIRON entry.

Cast iron material via MCIRON is supported for Small and Large Displacement Nonlinear Static Analysis and Nonlinear Transient Analysis. It is supported for CHEXA, CTETRA, CPENTA, and CPYRA elements (both first and second order).

In addition to regular result output, additional results from this material model are output if STRAIN(PLASTIC) I/O Option is specified. These results include, 6 plastic strain components, equivalent compressive plastic strain, equivalent tensile plastic strain, nodal results, and integration results are also available.

Heterogeneous materials such as composites, soil, solid foams, and polycrystals consist of one or more distinguishable constituents that have different mechanical properties. ¹² For example, it is known that fiber-reinforced composite materials have a stiffer and stronger phase known as the reinforcement and, a less stiff and weaker phase known as the matrix. The governing equations are well understood in one of the scales, either the micro-scale or the macro-scale. While general constitutive equations are applicable to solids under deformation in the macro-scale, continuum models are required in the micro-scale. In order to bridge the gap between the two scales, you need to consider models at different scales simultaneously. This can be achieved by multiscaling techniques.

Classical laminate theory can predict the overall properties of heterogeneous laminates; whereas, homogenization theory provides information on the individual phase properties and phase geometries. Knowledge of individual phases is important to characterize material damage and failure. Hence, multiscaling techniques open unique opportunities for modeling by focusing the material behavior on different levels of physical laws.

An overall flowchart for integration between the two software's is provided in Figure 7. Multiscale Designer has four pre-processing steps which are used to create material subroutines. These subroutines characterize the micro-scale material behavior in the linear and nonlinear domain. Later they are integrated with OptiStruct through the MATUSR material card for structural analysis.