Orthotropic damage model with fixed crack orientations for composites - CompDamMat

The model is designed for transversely isotropic elastic material defined by five elastic material constants. Typical example is a carbon fiber tow. Axis $1$ represents the material principal direction. The orthotropic material constants are defined as

$$\nu_{12} = \nu_{13},~\nu_{21}=\nu_{31},~\nu_{23}=\nu_{32},~E_{22}=E_{33}$$ (197)

$$G_{12} = G_{13}=G_{21}=G_{31},G_{23}=G_{32}$$ (198)

$$\frac{\nu_{12}=\nu_{13}}{E_{11}} = \frac{\nu_{21}=\nu_{31}}{E_{22}},~\frac{\nu_{31}=\nu_{21}}{E_{33}} = \frac{\nu_{13}=\nu_{12}}{E_{11}}$$ (199)

Material orientation on a finite element can be specified with lcs optional parameter. If unspecified, material orientation is the same as the global coordinate system. lcs array contains six numbers, where the first three numbers represent a directional vector of the local $x$-axis, and the next three numbers represent a directional vector of the local $y$-axis with the reference to the global coordinate system. The composite material is extended to 1D and is also suitable for beams and trusses. In such particular case, the lcs has no effect and the 1D element orientation is aligned with the global $xx$ component.

The index $p,~p\in\{11,22,33,23,31,12\}$ symbolizes six components of stress or strain vectors. The linear softening occurs after reaching a critical stress $f_{p,0}$, see Fig. 9. Orientation of cracks is assumed to be orthogonal and aligned with an orientation of material axes [3, pp.236]. The transverse isotropy is generally lost upon fracture, material becomes orthotropic and six damage parameters $d_p$ are introduced.

Figure 9: Implemented stress-strain evolution with damage for 1D case. Tension and compression are separated, but sharing the same damage parameter.

The compliance material matrix $\mathbf{H}$, in the secant form and including damage parameters, reads

$$\mathbf{H}= \left[ \begin{array}{cccccc} \frac{1}{(1-d_{11})E_{11... ...31})G_{31}} & 0\\ 0& 0&0 &0& 0& \frac{1}{(1-d_{12})G_{12}} \end{array} \right]$$ (200)

Damage occurs when any out of six stress tensor components exceeds a given strength $f_{p,0}$

$$\vert\sigma_p\vert \geq \vert f_{p,0}\vert$$ (201)

Positive and negative ultimate strengths can be generally different but share the same damage variable. At the point of damage initiation, see Fig. 9, one evaluates $\varepsilon_{p,E}$ and characteristic element length $l_p$, generally different for each damage mode. Given the fracture energy $G_{F,p}$, the maximum strain at zero stress $\varepsilon_{p,0}$ is computed

$$\varepsilon_{p,0} = \varepsilon_{p,E} + \frac{2G_{F,p}}{f_{p,0} l_p}$$ (202)

The point of damage initiation is never reached exactly, one needs to interpolate between the previous equilibrated step and current step to achieve objectivity.

The evolution of damage $d_p$ is based on the evolution of corresponding strain $\varepsilon_p$. A maximum achieved strain is stored in the variable $\kappa_p$. If $\varepsilon_p > \kappa_p$ the damage may grow so the corresponding damage variable $d_p$ may increase. Desired stress $\sigma'_p$ is evaluated from the actual strain $\varepsilon_p$

$$\sigma'_p = f_{p,0} \frac{\varepsilon_{p,0} - \varepsilon_p}{\varepsilon_{p,0}-\varepsilon_{p,E}}$$ (203)

and the calculation of damage variables $d_p$ stems from Eq. 200, for example

$$d_{11} = 1 - \frac{\sigma'_{11}}{E_{11}\left(\varepsilon_{11}+\frac{\nu_{21}}{E_{22}}\sigma_{22}+ \frac{\nu_{31}}{E_{33}}\sigma_{33} \right)}$$ (204)

$$d_{12} = 1 - \frac{\sigma'_{12}}{G_{12}\varepsilon_{12}}$$ (205)

Damage is always controled not to decrease. Fig. 10 shows a typical performance for this damage model in one direction.

The damage initiation is based on a trial stress. It becomes necessary for higher precision to skip a few first iteration, typically 5, and then to introduce damage. A parameter afterIter is designed for this purpose and MinIter forces a solver always to proceed certain amount of iterations.

allowSnapBack skips the checking of sufficient fracture energy for each direction. If not specified, all directions are checked to prevent snap-back which dissipates incorrect amount of energy.

Figure: Typical loading/unloading material performance for homogenized stress and strain in the direction `${2}$'. Note that one damage parameter is common for both tension and compression.

Table 45: Orthotropic damage model with fixed crack orientations for composites - summary.

Description	Orthotropic damage model with fixed crack orientations for composites
Record Format	CompDamMat num(in) # d(rn) # Exx(rn) # EyyEzz(rn) # nuxynuxz(rn) # nuyz(rn) # GxyGxz(rn) # Tension_f0_Gf(ra) # Compres_f0_Gf(ra) # [afterIter(in) #] [allowSnapBack(() #ia)]
Parameters	- num material model number
	- d material density
	- Exx Young's modulus for principal direction $xx$
	- EyyEzz Young's modulus in orthogonal directions to the principal direction $xx$
	- nuxynuxz Poisson's ratio in $xy$ and $xz$ directions
	- nuyz Poisson's ratio in $yz$ direction
	- GxyGxz shear modulus in $xy$ and $xz$ directions
	- Tension_f0_Gf array with six pairs of positive numbers. Each pair describes maximum stress in tension and fracture energy for each direction ($xx$, $yy$, $zz$, $yz$, $zx$, $xy$)
	- Compres_f0_Gf array with six pairs of numbers. Each pair describes maximum stress in compression (given as a negative number) and positive fracture energy for each direction ($xx$, $yy$, $zz$, $yz$, $zx$, $xy$)
Supported modes	3dMat, 1dMat
	- afterIter how many iterations must pass until damage is computed from strains, zero is default. User must ensure that the solver proceeds the minimum number of iterations.
	- allowSnapBack array to skip checking for snap-back. The array members are 1-6 for tension and 7-12 for compression components.