# NonlinearJ2

Nonlinear General J2 Plasticity Model

## Summary

This is an abstract material class thus cannot be used directly. This class defines a general plasticity model using J2 yielding criterion with associated flow rule and mixed hardening rule. The isotropic/kinematic hardening response can be customized.

To use this model, a derived class shall be defined first.

```cpp
class YourJ2 final : public NonlinearJ2 {
// class definition
}
```

The derived class only needs to implement four pure virtual methods that define the isotropic and kinematic hardening rules.

```cpp
virtual double compute_k(double) const = 0; // isotropic hardening
virtual double compute_dk(double) const = 0; // derivative isotropic
virtual double compute_h(double) const = 0; // kinematic hardening
virtual double compute_dh(double) const = 0; // derivative kinematic
```

All four methods take equivalent plastic strain as the input argument, on output, the corresponding quantities shall be provided.

The isotropic hardening function $$K(\bar\varepsilon\_p)$$ defines the isotropic hardening rule, there are some requirements:

1. $$K(\bar\varepsilon\_p)$$ should be non-negative,
2. $$K(\bar\varepsilon\_p=0)=\sigma\_y$$ where $$\sigma\_y$$ is the initial yielding stress.

There is no requirement for the kinematic hardening function $$H(\bar\varepsilon\_p)$$. Both hardening rules can coexist. However, to successfully solve the trial status, there is an additional constraint that shall be applied on the model:

$$
E+H'(\bar\varepsilon\_p)+K'(\bar\varepsilon\_p)\geqslant0\~\text{for all}\~\bar\varepsilon\_p
$$

Otherwise, the local Newton iteration will fail.

## Brief On Theory

The `NonlinearJ2` abstract class defines an associative plasticity framework using the von Mises yield criterion, which is defined as follows.

$$
F(\sigma,\bar\varepsilon\_p)=\sqrt{\dfrac{3}{2}(s-\beta(\bar\varepsilon\_p)):(s-\beta(\bar\varepsilon\_p))}-\sigma\_y(
\bar\varepsilon\_p)
$$

where $$\beta(\bar\varepsilon\_p)$$ is the back stress depends on the equivalent plastic strain $$\bar\varepsilon\_p$$ and $$\sigma\_y(\bar\varepsilon\_p)$$ is the yield stress. Note

$$
\sqrt{3J\_2}=\sqrt{\dfrac{3}{2}(s-\beta(\bar\varepsilon\_p)):(s-\beta(\bar\varepsilon\_p))}
$$

It is also called J2 plasticity model. A detailed discussion can be seen elsewhere. $$\beta(\bar\varepsilon\_p)=H( \bar\varepsilon\_p)$$ and $$\sigma\_y(\bar\varepsilon\_p)=K(\bar\varepsilon\_p)$$.

## History Layout

| location               | paramater                  |
| ---------------------- | -------------------------- |
| `initial_history(0)`   | accumulated plastic strain |
| `initial_history(1-6)` | back stress                |

## Kinematic Hardening

The back stress $$\beta(\bar{\varepsilon}\_p)$$ defines a kinematic hardening response. For example a linear kinematic hardening could be defined as:

$$
\beta(\bar\varepsilon\_p)=E\_K\bar\varepsilon\_p
$$

and the derivative

$$
\dfrac{\mathrm{d}\beta(\bar\varepsilon\_p)}{\mathrm{d}\bar\varepsilon\_p}=E\_K
$$

in which $$E\_K$$ is the kinematic hardening stiffness.

In this case, user shall override the corresponding two methods with such an implmentation.

```cpp
double SampleJ2::compute_h(const double p_strain) const { return e_kin * p_strain; }
double SampleJ2::compute_dh(const double) const { return e_kin; }
```

Of course, a nonlinear relationship could also be defined.

## Isotropic Hardening

The isotropic hardening is defined by function $$\sigma\_y(\bar\varepsilon\_p)$$. The value $$\sigma\_y(0)$$ should be the initial yield stress. Also, for a bilinear isotropic hardening response, user shall override the following two methods.

```cpp
double SampleJ2::compute_k(const double p_strain) const { return e_iso * p_strain + yield_stress; }
double SampleJ2::compute_dk(const double) const { return e_iso; }
```

Here another polynomial based isotropic hardening function is shown as an additional example. The function is defined as

$$
\sigma\_y=\sigma\_0(1+\sum\_{i=1}^{n}a\_i\bar\varepsilon\_p^i)
$$

where $$a\_i$$ are material constants. It is easy to see $$\sigma\_y|\_{\bar\varepsilon\_p=0}=\sigma\_0$$. The derivative is

$$
\dfrac{\mathrm{d}\sigma\_y}{\mathrm{d}\bar\varepsilon\_p}=\sigma\_0\sum\_{i=1}^{n}ia\_i\bar\varepsilon\_p^{i-1}
$$

To define such a hardening behavior, a vector shall be used to store all constants.

```cpp
// PolyJ2.h
const vec poly_para; // poly_para is a vector stores a_i
```

The methods could be written as follows.

```cpp
// PolyJ2.cpp
double PolyJ2::compute_k(const double p_strain) const {
 vec t_vec(poly_para.n_elem);

 t_vec(0) = 1.;
 for(uword I = 1; I < t_vec.n_elem; ++I) t_vec(I) = t_vec(I - 1) * p_strain;

 return yield_stress * dot(poly_para, t_vec);
}

double PolyJ2::compute_dk(const double p_strain) const {
 vec t_vec(poly_para.n_elem);

 t_vec(0) = 0.;
 t_vec(1) = 1.;
 for(uword I = 2; I < t_vec.n_elem; ++I) t_vec(I) = (double(I) + 1.) * pow(p_strain, double(I));

 return yield_stress * dot(poly_para, t_vec);
}
```

## Example

A few different models are shown as examples. User can define arbitrary models.

## Iso-error Map

The following example iso-error maps are obtained via the following script.

```py
from plugins import ErrorMap
# note: the dependency `ErrorMap` can be found in the following link
# https://github.com/TLCFEM/suanPan-manual/blob/dev/plugins/scripts/ErrorMap.py

young_modulus = 1e5
yield_stress = 100.0
hardening_ratio = 0.05

with ErrorMap(
    f"material BilinearJ2 1 {young_modulus} .2 {yield_stress} {hardening_ratio} 0.5",
    ref_strain=yield_stress / young_modulus,
    ref_stress=yield_stress,
    contour_samples=30,
) as error_map:
    error_map.contour("j2.uniaxial", center=(-4, 0), size=5, type={"rel", "abs"})
    error_map.contour("j2.biaxial", center=(-4, -4), size=5, type={"rel", "abs"})
```

![absolute error uniaxial](https://4006314410-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2F-MQ5rMqEBbzA9NgGETGX%2Fuploads%2Fgit-blob-700d2cd2e4f3db0117b99b9648c08c059d1a0d8a%2Fj2.uniaxial.abs.error.svg?alt=media) ![absolute error biaxial](https://4006314410-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2F-MQ5rMqEBbzA9NgGETGX%2Fuploads%2Fgit-blob-d2aea2836b4d77fb68fd1fc2e4cbb5f819ed3f68%2Fj2.biaxial.abs.error.svg?alt=media) ![relative error uniaxial](https://4006314410-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2F-MQ5rMqEBbzA9NgGETGX%2Fuploads%2Fgit-blob-0b5270bdf414bb88fd0c31fb1e3dfc19673baa39%2Fj2.uniaxial.rel.error.svg?alt=media) ![relative error biaxial](https://4006314410-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2F-MQ5rMqEBbzA9NgGETGX%2Fuploads%2Fgit-blob-14740da01162dd3a5fdb3accd3a3df85d328af83%2Fj2.biaxial.rel.error.svg?alt=media)