$\newcommand{\Dev}{\operatorname{Dev}}$ $\newcommand{\dev}{\operatorname{dev}}$ $\newcommand{\grad}{\operatorname{grad}}$ $\newcommand{\Grad}{\operatorname{Grad}}$ $\newcommand{\divg}{\operatorname{div}}$ $\newcommand{\Ei}{\operatorname{Ei}}$ $\newcommand{\tr}{\operatorname{tr}}$ $\newcommand{\Divg}{\operatorname{Div}}$ $\newcommand{\cay}{\operatorname{cay}}$ $\newcommand{\rot}{\operatorname{rot}}$

Curnier et al. [29] formulated a model for describing bimodular elastic solids exhibiting orthotropic material symmetry. This can be derived from the following hyperelastic strain-energy function: $$\begin{equation} W=\sum\limits _{a=1}^{3}\mu_{a}\mathbf{A}_{a}^{0}:\mathbf{E}^{2}+\frac{1}{2}\lambda_{aa}\left[\mathbf{A}_{a}^{0}:\mathbf{E}\right]\left(\mathbf{A}_{a}^{0}:\mathbf{E}\right)+\sum\limits _{b=1,\,b\ne a}^{3}\frac{1}{2}\lambda_{ab}\left(\mathbf{A}_{a}^{0}:\mathbf{E}\right)\left(\mathbf{A}_{b}^{0}:\mathbf{E}\right)\,,\label{eq434} \end{equation}$$ where $\mathbf{A}_{a}^{0}=\mathbf{a}_{a}^{0}\otimes\mathbf{a}_{a}^{0}$ is the structural tensor corresponding to one of the three mutually orthogonal planes of symmetry whose unit outward normal is $\mathbf{a}_{a}^{0}$ ($\mathbf{a}_{a}^{0}\cdot\mathbf{a}_{b}^{0}=\delta_{ab})$ . The bimodular response is described by $$\begin{equation} \lambda_{aa}\left[\mathbf{A}_{a}^{0}:\mathbf{E}\right]=\begin{cases} \lambda_{+aa} & \mathbf{A}_{a}^{0}:\mathbf{E}\geqslant0\\ \lambda_{-aa} & \mathbf{A}_{a}^{0}:\mathbf{E}<0 \end{cases}\,.\label{eq435} \end{equation}$$ The material constants are the three shear moduli $\mu_{a}$ , three tensile moduli $\lambda_{+aa}$ , three compressive moduli $\lambda_{-aa}$ , and three moduli $\lambda_{ab}$ ($b\ne a)$ , where $\lambda_{ba}=\lambda_{ab}$ . The second Piola-Kirchhoff stress can be derived from this strain energy density function: $$\begin{equation} \begin{aligned}\mathbf{S} & =\frac{\partial W}{\partial\mathbf{E}}=\sum\limits _{a=1}^{3}\mu_{a}\left(\mathbf{A}_{a}^{0}\cdot\mathbf{E}+\mathbf{E}\cdot\mathbf{A}_{a}^{0}\right)\\ & +\lambda_{aa}\left[\mathbf{A}_{a}^{0}:\mathbf{E}\right]\left(\mathbf{A}_{a}^{0}:\mathbf{E}\right)\mathbf{A}_{a}^{0}+\sum\limits _{b=1,\,b\ne a}^{3}\lambda_{ab}\left(\mathbf{A}_{a}^{0}:\mathbf{E}\right)\mathbf{A}_{b}^{0}\,. \end{aligned} \label{eq436} \end{equation}$$ The material elasticity tensor is then given by, $$\begin{equation} \begin{aligned} & \boldsymbol{\mathbb{C}}=\frac{\partial\mathbf{S}}{\partial\mathbf{E}}=\sum\limits _{a=1}^{3}\mu_{a}\left(\mathbf{A}_{a}^{0}\,\overline{\underline{\otimes}}\,\mathbf{I}+\mathbf{I}\,\overline{\underline{\otimes}}\,\mathbf{A}_{a}^{0}\right)\\ & +\lambda_{aa}\left[\mathbf{A}_{a}^{0}:\mathbf{E}\right]\mathbf{A}_{a}^{0}\otimes\mathbf{A}_{a}^{0}+\sum\limits _{b=1,\,b\ne a}^{3}\lambda_{ab}\mathbf{A}_{a}^{0}\otimes\mathbf{A}_{b}^{0} \end{aligned} \,.\label{eq437} \end{equation}$$ It is important to note that although this model is objective, it should only be used for small strains. For large strains, the response may be unrealistic. The Cauchy stress is $$\begin{equation} \begin{aligned}\boldsymbol{\sigma} & =J^{-1}\left(\sum\limits _{a=1}^{3}\frac{\mu_{a}}{2}\left(\mathbf{A}_{a}\cdot\left(\mathbf{b}-\mathbf{I}\right)+\left(\mathbf{b}-\mathbf{I}\right)\cdot\mathbf{A}_{a}\right)\right.\\ & \left.+\lambda_{aa}\left[K_{a}\right]K_{a}\mathbf{A}_{a}+\sum\limits _{b=1,\,b\ne a}^{3}\lambda_{ab}K_{a}\mathbf{A}_{b}\right)\,, \end{aligned} \label{eq438} \end{equation}$$ where $\mathbf{A}_{a}=\mathbf{F}\cdot\mathbf{A}_{a}^{0}\cdot\mathbf{F}^{T}$ and $K_{a}=\frac{1}{2}\left({\mathbf{A}_{a}:\mathbf{I}-1}\right)$ . The spatial elasticity tensor is $$\begin{equation} \boldsymbol{\mathcal{C}}=J^{-1}\left(\sum\limits _{a=1}^{3}\mu_{a}\left(\mathbf{A}_{a}\,\overline{\underline{\otimes}}\,\mathbf{b}+\mathbf{b}\,\overline{\underline{\otimes}}\,\mathbf{A}_{a}\right)+\lambda_{aa}\left[K_{a}\right]\mathbf{A}_{a}\otimes\mathbf{A}_{a}+\sum\limits _{b=1,\,b\ne a}^{3}\lambda_{ab}\mathbf{A}_{a}\otimes\mathbf{A}_{b}\right)\,.\label{eq439} \end{equation}$$ In the special case of cubic symmetry the number of material constants reduces to four, $$\begin{equation} \begin{aligned}\lambda_{+11} & =\lambda_{+22}=\lambda_{+33}\equiv\lambda_{+1}\\ \lambda_{-11} & =\lambda_{-22}=\lambda_{-33}\equiv\lambda_{-1}\\ \lambda_{12} & =\lambda_{23}=\lambda_{31}\equiv\lambda_{2}\\ \mu_{1} & =\mu_{2}=\mu_{3}\equiv\mu \end{aligned} \,.\label{eq440} \end{equation}$$