Wrinkling instabilities for biologically relevant fiber-reinforced composite materials with a case study of Neo-Hookean/Ogden–Gasser–Holzapfel bilayer

Abstract

Wrinkling is a ubiquitous surface phenomenon in many biological tissues and is believed to play an important role in arterial health. As arteries are highly nonlinear, anisotropic, multilayered composite systems, it is necessary to investigate wrinkling incorporating these material characteristics. Several studies have examined surface wrinkling mechanisms with nonlinear isotropic material relationships. Nevertheless, wrinkling associated with anisotropic constitutive models such as Ogden–Gasser–Holzapfel (OGH), which is suitable for soft biological tissues, and in particular arteries, still requires investigation. Here, the effects of OGH parameters such as fibers’ orientation, stiffness, and dispersion on the onset of wrinkling, wrinkle wavelength and amplitude are elucidated through analysis of a bilayer system composed of a thin, stiff neo-Hookean membrane and a soft OGH substrate subjected to compression. Critical contractile strain at which wrinkles occur is predicted using both finite element analysis and analytical linear perturbation approach. Results suggest that besides stiffness mismatch, anisotropic features associated with fiber stiffness and distribution might be used in natural layered systems to adjust wrinkling and subsequent folding behaviors. Further analysis of a bilayer system with fibers in the (x–y) plane subjected to compression in the x direction shows a complex dependence of wrinkling strain and wavelength on fiber angle, stiffness, and dispersion. This behavior is captured by an approximation utilizing the linearized anisotropic properties derived from OGH model. Such understanding of wrinkling in this artery wall-like system will help identify the role of wrinkling mechanisms in biological artery in addition to the design of its synthetic counterparts.

Appendices

Appendix 1: Perturbation analysis for wrinkling in neo-Hookean/Ogden–Gasser–Holzapfel Bilayer

Summary of notations

Symbols	Notation meaning
L	Length of the bilayer
H, h	Substrate and film thicknesses
$\lambda _x=\lambda _1, \lambda _2, \lambda _3$	Applied stretch in x, y, z directions, respectively
$\delta$	Perturbation amplitude
$\alpha$	Parameter determined in the eigenvalue analysis
$\lambda , k$	Wavelength and wave number: $k=\frac{2\pi }{\lambda }$
$\epsilon$	The applied strain: $\epsilon = \frac{\Delta L}{L}$
$\epsilon _{\mathrm{w}}$	Critical wrinkling strain
$k_1$, $k_2$	Fiber stiffness, fiber nonlinearity parameter
$\kappa$, $\theta$	Fiber dispersion, orientation
$L_i$	Unit vector of the direction of the $i\hbox {th}$ fiber family
$\sigma , P,S, N$	Cauchy, first and second Piola–Kirchhoff, nominal stresses
F, C, B	Deformation gradient, right and left Cauchy–Green tensor
$I_1, I_{4ii}=L_i^{\mathrm{T}}CL_i$	Invariants of the Cauchy–Green tensor
$\mu _{\mathrm{f}}, \mu _{\mathrm{s}}$	Shear modulus of film and substrate matrix
$u_{\mathrm{f}}, v_{\mathrm{f}}, u_{\mathrm{s}}, v_{\mathrm{s}}$	Displacement in the film and substrate

Bifurcation of a thin, stiff neo-Hookean layer on an OGH substrate

A bilayer composed of a thin, stiff incompressible neo-Hookean film attached to a soft, incompressible OGH substrate is subjected to compression as shown in Fig. 15.

Fig. 15

Thin stiff incompressible neo-Hookean film attached to a soft, incompressible OGH substrate with two symmetric families of fibers

Assume that the deformation is plane strain and uniform with the deformation state described in Eq. 8.

$$\begin{aligned} x=\lambda _x X \nonumber \\ y=\frac{1}{\lambda _x} Y \end{aligned}$$

(8)

where (X, Y) are coordinates in the undeformed configuration and (x, y) are the corresponding coordinates in the current configuration. $\lambda _x$ is the lateral stretch in the x direction.

Consider the perturbation of this homogeneous deformation state with the generic forms of perturbation shown in Eq. 9 for the case of incompressibility (Sun et al. 2012):

$$\begin{aligned} x& = \lambda _x X - \alpha \delta \lambda _x^2 \lambda \sin (k X) \hbox {e}^{\alpha k Y} \nonumber \\ y& = \frac{1}{\lambda _x} Y+\delta \lambda \cos (k X) \hbox {e}^{\alpha k Y} \nonumber \\ p& = p_0+\delta p_1 \cos (k X) \hbox {e}^{\alpha k Y} \end{aligned}$$

(9)

where $\delta<< 1$ is the perturbation amplitude parameter. $k=\frac{2\pi }{\lambda }$ is the wrinkle wave number, $\lambda$ is the undeformed-configurational wavelength of the perturbation. $\alpha$ is a parameter to be determined from the equilibrium condition. p is the hydrostatic pressure for incompressible solids. $p_0, p_1$ are to be determined from boundary conditions.

The deformation state of each layer in the bilayer due to the above perturbation is a linear combination of each layer ’s eigenmodes when each layer is considered separately and is subjected to the same perturbation. Therefore, in the followings, we first will consider the perturbation of each layer individually. A linear combination of the obtained eigenmodes will then be used to construct the deformation state of each layer in the bilayer system. Continuity conditions at the interface between the two layers and boundary conditions at the free surface are used next to construct an eigenvalue problem to determine the critical onset of wrinkling. Finally, the critical wrinkling strain is determined by numerically solving the resulting eigenvalue problem.

1.1 Analysis of eigenmodes of the neo-Hookean film

The analysis of eigenmodes of a neo-Hookean layer has been carried out in Sun et al. (2012), Cao and Hutchinson (2012), Stewart et al. (2016). Specifically, consider the neo-Hookean strain energy density function given in Eq. 10:

$$\begin{aligned} W_{\mathrm{NH}}=\mu _{\mathrm{f}} (I_1-3) \end{aligned}$$

(10)

First Piola–Kirchhoff stress can be determined from this strain energy function by taking the derivative with respect to the deformation gradients F and taking into account the hydrostatic pressure due to incompressibility, specifically:

$$\begin{aligned} P_{ij}=\frac{\partial W_{\mathrm{NH}}}{\partial F_{ij}} - p F_{ji}^{-1} \end{aligned}$$

(11)

Note that the deformation gradient can be determined from Eq. (9) for the perturbation state. Specifically,

$$\begin{aligned}&F_{11}=\frac{\partial {x}}{\partial {X}}=\lambda _x-\alpha \delta \lambda _x^2 \lambda k \left[ \cos (k X) \hbox {e}^{\alpha k Y}\right] \nonumber \\&F_{12}=\frac{\partial {x}}{\partial {Y}}=-\alpha ^2 \delta \lambda _x^2 \lambda k \left[ \sin (k X) \hbox {e}^{\alpha k Y}\right] \nonumber \\&F_{21}=\frac{\partial {y}}{\partial {X}}=-\delta \lambda k \left[ \sin (k X) \hbox {e}^{\alpha k Y}\right] \nonumber \\&F_{22}=\frac{\partial {y}}{\partial {Y}}=\frac{1}{\lambda _x}+\delta \lambda \alpha k \left[ \cos (k X) \hbox {e}^{\alpha k Y}\right] \nonumber \\&F_{13}=F_{31}=F_{23}=F_{32}=0, \quad F_{33}=1 \end{aligned}$$

(12)

Note that this prescribed deformation gradient already satisfies the incompressibility constraint $\hbox {det}(F)=1+O(\delta ^2)$.

The first invariant $I_1$ is the trace of the left Cauchy–Green tensor $B=F^{\mathrm{T}} F$, which is:

$$\begin{aligned} B& = \begin{bmatrix} F_{11} &{}\quad F_{12} &{}\quad 0 \\ F_{21} &{}\quad F_{22} &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 1 \end{bmatrix} \begin{bmatrix} F_{11} &{}\quad F_{21} &{}\quad 0 \\ F_{12} &{}\quad F_{22}&{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 1 \end{bmatrix} \nonumber \\& = \begin{bmatrix} F_{11}^2+F_{12}^2 &{}\quad F_{11}F_{21}+F_{12}F_{22} &{}\quad 0 \\ F_{21}F_{11}+F_{22}F_{12} &{}\quad F_{21}^2+F_{22}^2 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 1 \end{bmatrix} \end{aligned}$$

(13)

and

$$\begin{aligned} I_1=\hbox {trace}(B)=F_{11}^2+F_{12}^2+F_{21}^2+F_{22}^2+1 \end{aligned}$$

(14)

With the consideration of incompressibility $\hbox {det}(F)=1$ or $F_{11}F_{22}-F_{12}F_{21}=1$, the inverse of the deformation gradient $F^{-1}$ becomes:

$$\begin{aligned} F^{-1}& = \frac{1}{F_{11}F_{22}-F_{12}F_{21}} \begin{bmatrix} F_{22} &{}\quad -F_{12} &{}\quad 0 \\ -F_{21} &{}\quad F_{11} &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 1 \end{bmatrix} \nonumber \\& = \begin{bmatrix} F_{22} &{}\quad -F_{12} &{}\quad 0 \\ -F_{21} &{}\quad F_{11} &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 1 \end{bmatrix} \end{aligned}$$

(15)

Hence,

$$\begin{aligned} P_{11}& = \frac{\partial W_{\mathrm{NH}}}{\partial F_{11}} - p F^{-1}_{11}=\mu _{\mathrm{f}} \frac{\partial I_1}{\partial F_{11}} - p F^{-1}_{11} \nonumber \\& = 2\mu _{\mathrm{f}} F_{11}-\left[ p_0+\delta p_1 \cos (k X) \hbox {e}^{\alpha k Y} \right] F_{22} \nonumber \\ P_{12}& = \frac{\partial W_{\mathrm{NH}}}{\partial F_{12}} - p F^{-1}_{21}=\mu _{\mathrm{f}} \frac{\partial I_1}{\partial F_{12}} - p F^{-1}_{21} \nonumber \\& = 2\mu _{\mathrm{f}} F_{12}+\left[ p_0+\delta p_1 \cos (k X) \hbox {e}^{\alpha k Y} \right] F_{21} \nonumber \\ P_{21}& = \frac{\partial W_{\mathrm{NH}}}{\partial F_{21}} - p F^{-1}_{12}=\mu _{\mathrm{f}} \frac{\partial I_1}{\partial F_{21}} - p F^{-1}_{12} \nonumber \\& = 2\mu _{\mathrm{f}} F_{21}+\left[ p_0+\delta p_1 \cos (k X) \hbox {e}^{\alpha k Y} \right] F_{12} \nonumber \\ P_{22}& = \frac{\partial W_{\mathrm{NH}}}{\partial F_{22}} - p F^{-1}_{22}=\mu _{\mathrm{f}} \frac{\partial I_1}{\partial F_{22}} - p F^{-1}_{22} \nonumber \\& = 2\mu _{\mathrm{f}} F_{22}-\left[ p_0+\delta p_1 \cos (k X) \hbox {e}^{\alpha k Y} \right] F_{11} \end{aligned}$$

(16)

Note that for a more compact form, Eq. 16 can be written in the matrix form which can be easily derived from the second Piola Kirchhoff stress $S=2 \frac{\partial W_{\mathrm{NH}}}{\partial C}$. Specifically,

$$\begin{aligned}&S=2 \mu _{\mathrm{f}} \frac{\partial I_{1}}{\partial C}=2 \mu _{\mathrm{f}} I \nonumber \\&P=FS - p F^{-T}=2 \mu _{\mathrm{f}} F - p F^{-T} \end{aligned}$$

(17)

which are consistent with the stresses using index notations in Eq. 16.

Substituting Eq. 12 into Eq. 16, the calculations result in the following formulae for first Piola–Kirchhoff stress:

$$\begin{aligned} P_{11}& = \left( 2\mu _{\mathrm{f}} \lambda _x - p_0\frac{1}{\lambda _x}\right) \nonumber \\&\quad -\, \delta \left\{ 2\mu _{\mathrm{f}} k\alpha \lambda _x^2 \lambda +p_0 k \alpha \lambda +\frac{p_1}{\lambda _x}\right\} \cos (kX) \hbox {e}^{\alpha k Y} \nonumber \\&\quad +\, O(\delta ^2) \nonumber \\ P_{12}& = -\delta \left( 2\mu _{\mathrm{f}} k \alpha ^2 \lambda _x^2 \lambda + p_0 k \lambda \right) \sin (kX) \hbox {e}^{\alpha k Y} + O(\delta ^2) \nonumber \\ P_{21}& = -\delta \left( 2\mu _{\mathrm{f}} k \lambda + p_0 k \alpha ^2 \lambda _x^2 \lambda \right) \sin (kX)\hbox {e}^{\alpha k Y}+ O(\delta ^2) \nonumber \\ P_{22}& = \left( \frac{2\mu _{\mathrm{f}}}{\lambda _x}-p_0\lambda _x\right) \nonumber \\&\quad +\,\delta \left( 2\mu _{\mathrm{f}} k \alpha \lambda + p_0 k \alpha \lambda _x^2 \lambda - \lambda _x p_1\right) \cos (kX) \hbox {e}^{\alpha k Y} \nonumber \\&\quad +\, O(\delta ^2) \end{aligned}$$

(18)

The zeroth order ($\delta =0$) solution is obtained from Eq. 18, and from the boundary condition for the single layer $P_{22}^{0}=0$, it is shown that $p_0=\frac{2\mu _{\mathrm{f}}}{\lambda _x^2}$

Substituting the stresses into the two following equilibrium equations, the following 2 equations are obtained:

$$\begin{aligned} P_{11,X}+P_{12,Y}& = \left\{ 4\pi \mu _{\mathrm{f}} \alpha \lambda _x^2 \left( 1-\alpha ^2\right) +\frac{p_1}{\lambda _x}\right\} \nonumber \\&\quad \delta k \sin (kX) \hbox {e}^{\alpha k Y} + O(\delta ^2)=0 \nonumber \\ P_{12,X}+P_{22,Y}& = \left( -4\pi \mu _{\mathrm{f}} + 4 \pi \mu _{\mathrm{f}} \alpha ^2 - \lambda _x \alpha p_1\right) \nonumber \\&\quad \delta k \cos (kX) \hbox {e}^{\alpha k Y} + O(\delta ^2) = 0 \end{aligned}$$

(19)

which leads to:

$$\begin{aligned}&4\pi \mu _{\mathrm{f}} \alpha \lambda _x^2 \left( 1-\alpha ^2\right) +\frac{p_1}{\lambda _x}=0 \nonumber \\&-\,4\pi \mu _{\mathrm{f}} + 4 \pi \mu _{\mathrm{f}} \alpha ^2 - \lambda _x \alpha p_1= 0 \end{aligned}$$

(20)

From the first part of Eq. 20, $p_1=-4\pi \mu _{\mathrm{f}} \alpha \lambda _x^3 (1-\alpha ^2)$. Substituting this $p_1$ into the second part of Eq. 20 yields a fourth-order equation of $\alpha$: $4\pi \mu _{\mathrm{f}} (\lambda _x^4\alpha ^2-1)(1-\alpha ^2)=0$. Solving this equation gives four solutions of $\alpha$ and hence 4 pairs of solutions $\left( \alpha _i, p_i\right)$ corresponding with 4 eigenvalues: $\alpha =1, \alpha =-1, \alpha =1/\lambda _x^2, \alpha =-1/\lambda _x^2$. Substituting each pair of solutions into Eqs. (9) and (18) provides an eigenmode and its stress state for the single neo-Hookean layer subjected to perturbation. Specifically, from Eq. (9), $\left( u_{fi},v_{fi}\right) , i=\overline{1,4}$ are obtained for the deformation where $u=x-\lambda _x X, v=y-Y/\lambda _x$. From Eq. (18), tangential and normal tractions $N_{fi_{21}}, N_{fi_{22}}, i=\overline{1,4}$ are obtained from the nominal stresses, respectively. Nominal stress is determined as $N= P^{\mathrm{T}}$.

1.2 Analysis of eigenmodes of the OGH substrate

Consider the OGH substrate with the strain energy density function given in Eq. 21.

$$\begin{aligned} W_{\mathrm{OGH}}& = \mu _{\mathrm{s}} (I_1-3)+\frac{k_1}{2k_2}\sum _{m=1}^{N}\left\{ \hbox {exp}\left[ k_2 E_{\mathrm{m}}^2\right] -1\right\} \nonumber \\ E_{\mathrm{m}}& = \kappa (I_1-3)+(1-3\kappa )(I_\mathrm{4mm}-1) \end{aligned}$$

(21)

where the invariant $I_\mathrm{4mm}=L_{\mathrm{m}}.(C L_{\mathrm{m}})$, $L_{\mathrm{m}}$ is the unit direction of the fiber family m th, $C=F^{\mathrm{T}}F$ is the right Cauchy Green tensor.

With the same approach as in Sect. 1.1, the first Piola–Kirchhoff stress can be determined by taking the derivative of the energy with respect to the deformation gradient.

$$\begin{aligned} P_{ij}=\frac{\partial W_{\mathrm{OGH}}}{\partial F_{ij}} - p F_{ji}^{-1} \end{aligned}$$

(22)

or from the second Piola–Kirchhoff stress S:

$$\begin{aligned} P=FS- p F^{-T}=2F\frac{\partial W_{\mathrm{OGH}}}{\partial C} - p F^{-T} \end{aligned}$$

(23)

where the second Piola–Kirchhoff stress is:

$$\begin{aligned} S& = 2\frac{\partial W_{\mathrm{OGH}}}{\partial C}=2\mu _{\mathrm{s}}\frac{\partial {I_1}}{\partial {C}} \nonumber \\&\quad +\frac{k_1}{k_2}\sum _{m=1}^2 \left\{ 2k_2E_{\mathrm{m}} \hbox {exp}\left[ k_2E_{\mathrm{m}}^2\right] \frac{\partial {E_i}}{\partial {C}}\right\} \nonumber \\ S& = 2\mu _{\mathrm{s}} I+2k_1\sum _{m=1}^2 E_{\mathrm{m}} \hbox {exp}\left[ k_2 E_{\mathrm{m}}^2\right] \nonumber \\&\quad \left\{ \kappa \frac{\partial {I_1}}{\partial {C}}+(1-3\kappa ) \frac{\partial {I_{4mm}}}{\partial {C}}\right\} \nonumber \\ S& = 2\mu _{\mathrm{s}} I+ \sum _{m=1}^2 2 k_1 \kappa E_{\mathrm{m}} \hbox {exp}\left[ k_2 E_{\mathrm{m}}^2\right] I \nonumber \\&\quad + \sum _{m=1}^2 2 k_1 (1-3\kappa ) E_{\mathrm{m}} \hbox {exp}\left[ k_2 E_{\mathrm{m}}^2\right] L_{\mathrm{m}} \otimes L_{\mathrm{m}} \end{aligned}$$

(24)

Thus, the first Piola–Kirchhoff stress becomes:

$$\begin{aligned} P& = -PF^{-T}+2\mu _{\mathrm{s}} F+ \sum _{m=1}^2 2 k_1 \kappa E_{\mathrm{m}} \hbox {exp}\left[ k_2 E_{\mathrm{m}}^2\right] F \nonumber \\&+\sum _{m=1}^2 2 k_1 (1-3\kappa ) E_{\mathrm{m}} \hbox {exp}\left[ k_2 E_{\mathrm{m}}^2\right] F\left( L_{\mathrm{m}} \otimes L_{\mathrm{m}}\right) \end{aligned}$$

(25)

With the given deformation gradient in Eq. (9), components of the first Piola–Kirchhoff stress corresponding to the OGH layer are specified in Eq. 25. Hence, from the zeroth-order solution with the boundary condition $P_{22}^{0}=0$, $p_0$ is determined.

Substituting these stresses into the two equilibrium equations (in Eq. 19) and applying the same solution method as in Sect. 1.1, a fourth-order equation in terms of $\alpha$ is again obtained. Hence, 4 pairs of solutions $\left( \alpha _i, p_i\right)$ are determined. However, for the substrate, as the undulation dies down as $Y \rightarrow -\infty$, only 2 solutions with positive values of $\alpha$ are used to construct eigenmodes for the substrate. Specifically, from Eq. (9), $\left( u_{si},v_{si}\right) , i=\overline{1,2}$ are obtained for the deformation where $u=x-\lambda _x X, v=y-Y/\lambda _x$. From Eq. (25), tangential and normal tractions $N_{si_{21}}, N_{si_{22}}, i=\overline{1,2}$ are obtained from the nominal stresses, respectively. Nominal stress is determined as $N= P^{\mathrm{T}}$.

Note that for a neo-Hookean layer, the 4 eigenvalues $\alpha$ are all real values. However, solving the fourth-order equation of eigenvalues $\alpha$ for the OGH substrate can result in complex solutions. If the eigenvalues $\alpha$ for the OGH substrate are complex, the four eigenvalues will correspond to two pairs of complex conjugate. As the undulation must vanish when $Y \rightarrow -\infty$, the pair of complex conjugate with the positive real part is chosen for constructing the eigenmodes for the substrate. Actually, as this is a pair of conjugate eigenvalues $\alpha$, only one is needed. With this eigenvalue, it is straightforward to compute $\left( u_{\mathrm{si}},v_{\mathrm{si}}\right)$ and $N_{{{\text{si}}_{{21}} }} ,N_{{{\text{si}}_{{22}} }}$ which are the corresponding deformation and stresses. As $\alpha$ is a complex value, these deformation and stress fields are also complex. The real parts and the complex parts of the fields are used now to construct the 2 eigenmodes of the OGH substrate.

Due to the cumbersome formulae, all the calculations for the stresses, equilibrium equations, and eigenvalue problems are implemented in Matlab.

1.3 Linear combination of eigenmodes for the bilayer system

Recall: $\left( u_{fi},v_{fi}\right) , i=\overline{1,4}$ and $N_{fi_{21}}, N_{fi_{22}}, i=\overline{1,4}$ as the deformation and stresses, respectively, associated with 4 eigenmodes of the neo-Hookean layer. $\left( u_{si},v_{si}\right) , i=\overline{1,2}$ and $N_{si_{21}}, N_{si_{22}}, i=\overline{1,2}$ as the deformation and nominal stresses, respectively, associated with 2 eigenmodes of the OGH layer.

For a bilayer composed of a neo-Hookean film attached to an OGH substrate subjected to the perturbation in Eq. 9, the deformation in each layer is a linear combination of its eigenmodes. In other words, the deformation in the neo-Hookean film can be written as:

$$\begin{aligned} u_{\mathrm{f}} &= A_1 u_{\mathrm{f1}}+A_2 u_{\mathrm{f2}}+A_3 u_{\mathrm{f3}}+A_4 u_{\mathrm{f4}} \nonumber \\ v_{\mathrm{f}} &= A_1 v_{\mathrm{f1}}+A_2 v_{\mathrm{f2}}+A_3 v_{\mathrm{f3}}+A_4 v_{\mathrm{f4}} \nonumber \\ N_{{\mathrm{f}}_{21}} &= A_1 N_{{\mathrm{f1}}_{21}}+A_2 N_{{\mathrm{f2}}_{21}}+A_3 N_{{\mathrm{f3}}_{21}}+A_4 N_{{\mathrm{f4}}_{21}} \nonumber \\ N_{{\mathrm{f}}_{22}} &= A_1 N_{{\mathrm{f1}}_{22}}+A_2 N_{{\mathrm{f2}}_{22}}+A_3 N_{{\mathrm{f3}}_{22}}+A_4 N_{{\mathrm{f4}}_{22}} \end{aligned}$$

(26)

The deformation in the OGH substrate can be written as:

$$\begin{aligned} u_{\mathrm{s}}& = A_5 u_{\mathrm{s1}}+A_6 u_{\mathrm{s2}} \nonumber \\ v_{\mathrm{s}}& = A_5 v_{\mathrm{s1}}+A_6 v_{\mathrm{s2}} \nonumber \\ N_{{\mathrm{s}}_{21}}& = A_5 N_{{\mathrm{s1}}_{21}}+A_6 N_{{\mathrm{s2}}_{21}} \nonumber \\ N_{{\mathrm{s}}_{22}}& = A_5 N_{{\mathrm{s1}}_{22}}+A_6 N_{{\mathrm{s2}}_{22}} \end{aligned}$$

(27)

where $A_1,A_2,A_3,A_4,A_5,A_6$ are constant parameters.

1.4 Continuity and boundary conditions—eigenvalue problem for critical strain

At the interface between the film and the substrate, $Y=0$, continuity in displacements and tractions are enforced which can be written as follows:

$$\begin{aligned} u_{\mathrm{f}}(Y=0)& = u_{\mathrm{s}}(Y=0) \nonumber \\ v_{\mathrm{f}}(Y=0)& = v_{\mathrm{s}}(Y=0) \nonumber \\ N_{{\mathrm{f}}_{21}}(Y=0)& = N_{{\mathrm{s}}_{21}}(Y=0) \nonumber \\ N_{{\mathrm{f}}_{22}}(Y=0)& = N_{{\mathrm{s}}_{22}}(Y=0) \end{aligned}$$

(28)

At the stress-free face $Y=h$ of the neo-Hookean film, two boundary conditions for tractions are obtained:

$$\begin{aligned} N_{{\mathrm{f}}_{21}}(Y=h)& = 0\nonumber \\ N_{{\mathrm{f}}_{22}}(Y=h)& = 0 \end{aligned}$$

(29)

1.5 Critical wrinkling strain determination

By substituting the linear combinations in Eqs. 26, 27 into the system of continuity and boundary conditions in Eqs. 28 and 29, a system of the following forms is obtained:

$$\begin{aligned} F\left( \lambda _x,kh\right) \begin{bmatrix} A_1\\ A_2\\ A_3\\ A_4\\ A_5\\ A_6 \end{bmatrix} =0 \end{aligned}$$

(30)

The critical wrinkling strain, which is minimized over all kh is determined from the nonlinear equations $\det (F)=0$. This is solved numerically in Matlab.

Appendix 2: Linearization of OGH model for material properties

Figure 16 plots the prediction of $\frac{{(kh)}^2}{\epsilon _{\mathrm{w}}}$ ratio for bilayers of a neo-Hookean film bonded to an OGH substrate.

Fig. 16

Ratio$\frac{{(kh)}^2}{\epsilon _{\mathrm{w}}}$ versus fiber angle $\theta$ for the case of fiber dispersion $\kappa =0$ in (x–y) fiber plane

When the fiber stiffness $k_1=0$, the substrate becomes neo-Hookean. The ratio approaches the value of 4 for $\mu _{\mathrm{f}}>> \mu _{\mathrm{s}}$, which is demonstrated in Fig. 16 for the mismatch modulus ratio $\mu _{\mathrm{f}}/\mu _{\mathrm{s}}=1000$. At lower mismatch ratio, the value of $\frac{{(kh)}^2}{\epsilon _{\mathrm{w}}}$ is slightly less than 4. Specifically, with a mismatch modulus ratio $\mu _{\mathrm{f}}/\mu _{\mathrm{s}}$ of 100 as considered in this paper, $\frac{{(kh)}^2}{\epsilon _{\mathrm{w}}} \sim 3.7$. When fiber stiffness is nonzero, the renormalization of the substrate stiffness must change the effective value of the substrate stiffness and lead to deviations of this ratio from the constant value. At low fiber stiffness $k_1/\mu _{\mathrm{s}}=2, k_1/\mu _{\mathrm{s}}=5$, the ratio remains constant around the value of the corresponding neo-Hookean bilayer with the same modulus mismatch $\mu _{\mathrm{f}}/\mu _{\mathrm{s}}=100$. For higher fiber stiffness $k_1/\mu _{\mathrm{s}}=10$, the ratio shows some deviations from this constant trend, especially at high values of fiber angle $\theta$. Note that as angle $\theta$ increases, the transverse direction y also has higher stiffness. The deviation, therefore, might be attributed to the effect of orthotropic material properties associated with OGH model which become more significant as fibers are stiffer and oriented in the transverse (y) direction. These properties are derived through linearization as follows.

Ogden–Gasser–Holzapfel substrate with strain energy density function:

$$\begin{aligned} W_{\mathrm{OGH}}& = \mu (I_1-3)+\frac{k_1}{2k_2}\sum _{i=1}^{N}\left\{ \hbox {exp}\left[ k_2 E_i^2\right] -1\right\} \nonumber \\ E_i& = \kappa (I_1-3)+(1-3\kappa )(I_{4ii}-1) \nonumber \\ I_{4ii}& = L_i^{\mathrm{T}} C L_i \end{aligned}$$

(31)

where $C=F^{\mathrm{T}}F$ is the right Cauchy–Green tensor, F is the deformation gradient, $L_i$ is the unit vector of the orientation of the $i\hbox {th}$ fiber family. Here, $N=2$ corresponds to two fiber families.

The second Piola–Kirchhoff stress:

$$\begin{aligned} S& = 2\frac{\partial W_{\mathrm{OGH}}}{\partial C} \nonumber \\& = 2\mu \frac{\partial {I_1}}{\partial {C}}+\frac{k_1}{k_2}\sum _{i=1}^2 \left\{ 2k_2E_i \hbox {exp}\left[ k_2E_i^2\right] \frac{\partial {E_i}}{\partial {C}}\right\} \nonumber \\ S& = 2\mu I+2k_1\sum _{i=1}^2 E_i \hbox {exp}\left[ k_2 E_i^2\right] \nonumber \\&\quad \left\{ \kappa \frac{\partial {I_1}}{\partial {C}}+(1-3\kappa ) \frac{\partial {I_{4i}}}{\partial {C}}\right\} \nonumber \\ S& = 2\mu I+ \sum _{i=1}^2 2 k_1 \kappa E_i \hbox {exp}\left[ k_2 E_i^2\right] I \nonumber \\&\quad + \sum _{i=1}^2 2 k_1 (1-3\kappa ) E_i \hbox {exp}\left[ k_2 E_i^2\right] L_i \otimes L_i \end{aligned}$$

(32)

Cauchy stress for incompressible case:

$$\begin{aligned} \sigma& = FSF^{\mathrm{T}}=-pI+2\mu B+ \sum _{i=1}^2 2 k_1 \kappa E_i \hbox {exp}\left[ k_2 E_i^2\right] B \nonumber \\&\quad + \sum _{i=1}^2 2 k_1 (1-3\kappa ) E_i \hbox {exp}\left[ k_2 E_i^2\right] FL_i \otimes FL_i \end{aligned}$$

(33)

where $B=FF^{\mathrm{T}}$ is the left Cauchy Green tensor.

For two family fibers lying in x–y plane: $L_1=[c, s, 0]^{\mathrm{T}}, L_2=[c, -s, 0]^{\mathrm{T}}, c=\cos (\theta ), s=\sin (\theta )$ where $\theta$ is the fiber angle with respect to x-axis.

1.1 Determine longitudinal moduli and Poisson’s ratios

Consider a block made of OGH material being subjected to tension in one direction and free to expand in the other two directions. The deformation gradient F, left Cauchy–Green tensor B, right Cauchy–Green tensor C are described as follows:

$$\begin{aligned} F=\begin{bmatrix} \lambda _1 &{}\quad 0 &{}\quad 0 \\ 0 &{}\quad \lambda _2 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad \lambda _3 \end{bmatrix}, \quad B=C=\begin{bmatrix} \lambda _1^2 &{}\quad 0 &{}\quad 0 \\ 0 &{}\quad \lambda _2^2 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad \lambda _3^2 \end{bmatrix} \end{aligned}$$

(34)

where the three stretch ratios are related by incompressibility restriction: $\lambda _1\lambda _2\lambda _3=1$. Other quantities in Eq. (33) for computing Cauchy stresses become:

$$\begin{aligned} FL_1& = \begin{bmatrix} \lambda _1 &{}\quad 0 &{}\quad 0 \\ 0 &{}\quad \lambda _2 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad \lambda _3 \end{bmatrix} \begin{bmatrix} c\\ s\\ 0 \end{bmatrix} = \begin{bmatrix} c\lambda _1\\ s\lambda _2\\ 0 \end{bmatrix},\nonumber \\ FL_2& = \begin{bmatrix} \lambda _1 &{}\quad 0 &{}\quad 0 \\ 0 &{}\quad \lambda _2 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad \lambda _3 \end{bmatrix} \begin{bmatrix} c\\ -s\\ 0 \end{bmatrix} = \begin{bmatrix} c\lambda _1\\ -s\lambda _2\\ 0 \end{bmatrix} \end{aligned}$$

(35)

$$\begin{aligned} I_{41}& = \begin{bmatrix} c&\quad s&\quad 0 \end{bmatrix} \begin{bmatrix} \lambda _1^2 &{}\quad 0 &{}\quad 0 \\ 0 &{}\quad \lambda _2^2 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad \lambda _3^2 \end{bmatrix} \begin{bmatrix} c\\ s\\ 0 \end{bmatrix} \nonumber \\& = c^2\lambda _1^2+s^2\lambda _2^2 \end{aligned}$$

(36)

$$\begin{aligned} I_{42}& = \begin{bmatrix} c&\quad -s&\quad 0 \end{bmatrix} \begin{bmatrix} \lambda _1^2 &{}\quad 0 &{}\quad 0 \\ 0 &{}\quad \lambda _2^2 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad \lambda _3^2 \end{bmatrix} \begin{bmatrix} c\\ -s\\ 0 \end{bmatrix} \nonumber \\& = c^2\lambda _1^2+s^2\lambda _2^2 \end{aligned}$$

(37)

$$\begin{aligned} E_1& = E_2=\kappa (I_1-3) \nonumber \\&\quad + (1-3\kappa )(c^2\lambda _1^2+s^2\lambda _2^2-1) \end{aligned}$$

(38)

$$\begin{aligned} FL_1 \otimes FL_1& = \begin{bmatrix} c\lambda _1 \\ s\lambda _2 \\ 0 \end{bmatrix} \begin{bmatrix} c\lambda _1&\quad s\lambda _2&\quad 0 \end{bmatrix} \nonumber \\& = \begin{bmatrix} c^2\lambda _1^2 &{}\quad cs\lambda _1\lambda _2 &{}\quad 0\\ cs\lambda _1\lambda _2 &{}\quad s^2\lambda _2^2 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 0 \end{bmatrix} \end{aligned}$$

(39)

$$\begin{aligned} FL_2 \otimes FL_2& = \begin{bmatrix} c\lambda _1 \\ - s\lambda _2 \\ 0 \end{bmatrix} \begin{bmatrix} c\lambda _1&\quad -s\lambda _2&\quad 0 \end{bmatrix} \nonumber \\& = \begin{bmatrix} c^2\lambda _1^2 &{}\quad -cs\lambda _1\lambda _2 &{}\quad 0\\ -cs\lambda _1\lambda _2 &{}\quad s^2\lambda _2^2 &{}\quad 0 \\ 0 &{} 0 &{} 0 \end{bmatrix} \end{aligned}$$

(40)

Substituting Eqs. (34–40) into Eq. (33), three components of Cauchy stresses for the block being pulled in one direction:

$$\begin{aligned} \sigma _{11}& = -p+2\mu \lambda _1^2+4 \kappa k_1 E_1 \hbox {exp}\left[ k_2 E_1^2\right] \lambda _1^2 \nonumber \\&\quad +\, 4 (1-3\kappa ) k_1 E_1 \hbox {exp}\left[ k_2 E_1^2\right] c^2 \lambda _1^2 \nonumber \\ \sigma _{22}& = -p+2\mu \lambda _2^2+4 \kappa k_1 E_1 \hbox {exp}\left[ k_2 E_1^2\right] \lambda _2^2 \nonumber \\&\quad +\, 4 (1-3\kappa ) k_1 E_1 \hbox {exp}\left[ k_2 E_1^2\right] s^2 \lambda _2^2 \nonumber \\ \sigma _{33}& = -p+2\mu \lambda _3^2+4 \kappa k_1 E_1 \hbox {exp}\left[ k_2 E_1^2\right] \lambda _3^2 \end{aligned}$$

(41)

1.1.1 Stretching along the x direction to determine $E_x, \nu _{xy}, \nu _{xz}$

Note that the block is pulled in x direction and is free to expand in y and z directions, therefore:

$$\begin{aligned} \sigma _{11} \ne 0, \sigma _{22}=0, \sigma _{33}=0 \end{aligned}$$

(42)

thus, using $\lambda _3=\frac{1}{\lambda _1\lambda _2}$, they can be written:

$$\begin{aligned} \sigma _{11}& = \sigma _{11}-\sigma _{33}=2\mu \left\{ \lambda _1^2- \frac{1}{\lambda _1^2\lambda _2^2}\right\} \nonumber \\&\quad +\,4\kappa k_1 E_1 \hbox {exp}[k_2 E_1^2] \left[ \lambda _1^2- \frac{1}{\lambda _1^2\lambda _2^2}\right] \nonumber \\&\quad +\,4(1-3\kappa )k_1E_1 \hbox {exp}[k_2 E_1^2]c^2\lambda _1^2 \nonumber \\ \sigma _{22}& = \sigma _{22}-\sigma _{33}=2\mu \left\{ \lambda _2^2- \frac{1}{\lambda _1^2\lambda _2^2}\right\} \nonumber \\&\quad +\,4\kappa k_1 E_1 \hbox {exp}[k_2 E_1^2] \left[ \lambda _2^2-\frac{1}{\lambda _1^2\lambda _2^2}\right] \nonumber \\&\quad +\,4(1-3\kappa )k_1E_1 \hbox {exp}[k_2 E_1^2]s^2\lambda _2^2 \end{aligned}$$

(43)

Consider small deformation regime, the strains are small and their high order terms can be neglected. Thus, the following approximations can be used to approximate the stresses:

$$\begin{aligned}&\lambda _1=1+\epsilon _1, \lambda _2=1+\epsilon _2, \lambda _1^2=1+2\epsilon _1, \lambda _2^2=1+2\epsilon _2 \nonumber \\&\frac{1}{\lambda _1^2}=1-2\epsilon _1, \frac{1}{\lambda _2^2}=1-2\epsilon _2, \frac{1}{\lambda _1^2\lambda _2^2}=1-2\epsilon _1-2\epsilon _2 \nonumber \\&I_1=\lambda _1^2+\lambda _2^2+\lambda _3^2=3 \nonumber \\&I_{41}=I_{42}=c^2\lambda _1^2+s^2\lambda _2^2=1+2c^2\epsilon _1+2s^2\epsilon _2 \nonumber \\&E_1=E_2=(1-3\kappa )\left( 2c^2\epsilon _1+2s^2\epsilon _2\right) \nonumber \\&\hbox {exp}\left[ k_2 E_1^2\right] =1 \end{aligned}$$

(44)

Substituting the approximations in Eq. (44) into Eq. (43), the stresses become:

$$\begin{aligned} \sigma _{11}& = 2\mu (4\epsilon _1+2\epsilon _2)+8k_1(1-3\kappa )^2 (c^4\epsilon _1+c^2s^2\epsilon _2) \nonumber \\ \sigma _{22}& = 2\mu (4\epsilon _2+2\epsilon _1)+8k_1(1-3\kappa )^2 (c^2s^2\epsilon _1+s^4\epsilon _2) \end{aligned}$$

(45)

since $\sigma _{22}=0$, so we have:

$$\begin{aligned} \frac{\epsilon _2}{\epsilon _1}=-\nu _{xy}=-\frac{4\mu +8k_1(1-3\kappa )^2 c^2 s^2}{8\mu +8k_1(1-3\kappa )^2s^4} \end{aligned}$$

(46)

Note that we are pulling in the x direction, so this ratio between the two strains gives the Poisson’s ratio $\nu _{xy}$. It can be seen that for the case of isotropic material, i.e., $k_1=0$ or $\kappa =1/3$, this ratio is equal to 0.5 which is the Poisson’s ratio of isotropic incompressible material.

By substituting $\epsilon _2$ in terms of $\epsilon _1$ into $\sigma _{11}$, we have:

$$\begin{aligned} \sigma _{11}=\frac{6\mu ^2+8k_1 \mu (1-3\kappa )^2(1-3c^2s^2)}{\mu +k_1(1-3\kappa )^2s^4} \epsilon _1 \end{aligned}$$

(47)

The longitudinal modulus in the x-direction, thus, can be obtained:

$$\begin{aligned} E_x=\frac{\sigma _{11}}{\epsilon _{11}}=\frac{6\mu ^2+8k_1 \mu (1-3\kappa )^2(1-3c^2s^2)}{\mu +k_1(1-3\kappa )^2s^4} \end{aligned}$$

(48)

Again, for $k_1=0$ or $\kappa =1/3$, this modulus reduces to $E_x=6\mu$ which is the Young modulus value of isotropic, incompressible material.

Note that:

$$\begin{aligned}&\lambda _3=\frac{1}{\lambda _1\lambda _2}=\frac{1}{(1+\epsilon _1)(1+\epsilon _2)} =1-\epsilon _1-\epsilon _2=1+\epsilon _3\nonumber \\&\epsilon _1+\epsilon _2=-\epsilon _3 \implies 1+\frac{\epsilon _2}{\epsilon _1} =-\frac{\epsilon _3}{\epsilon _1}\nonumber \\&\implies \nu _{xz}=-\frac{\epsilon _3}{\epsilon _1}=\frac{4\mu +8k_1(1-3\kappa )^2 s^2(s^2-c^2)}{8\mu +8k_1(1-3\kappa )^2s^4} \end{aligned}$$

(49)

1.1.2 Stretching along the y direction to determine $E_y, \nu _{yx}, \nu _{yz}$

With the same approach, the modulus $E_y$ in y-direction and Poisson’s ratio $\nu _{yx}$ can be derived by subjecting the block to the tension in y-direction. Specifically,

$$\begin{aligned} \sigma _{11} = 0, \sigma _{22} \ne 0, \sigma _{33}=0 \end{aligned}$$

(50)

Thus:

$$\begin{aligned} \sigma _{11}& = \sigma _{11}-\sigma _{33}=\left\{ 8\mu +8k_1(1-3\kappa )^2 c^4\right\} \epsilon _1 \nonumber \\&\quad +\,\left\{ 4\mu +8k_1(1-3\kappa )^2 c^2s^2\right\} \epsilon _2=0 \nonumber \\ \sigma _{22}& = \sigma _{22}-\sigma _{33}=\left\{ 4\mu +8k_1(1-3\kappa )^2 c^2s^2\right\} \epsilon _1 \nonumber \\&\quad +\,\left\{ 8\mu +8k_1(1-3\kappa )^2s^4\right\} \epsilon _2 \end{aligned}$$

(51)

The Poisson’s ratio $\nu _{yx}$ is therefore,

$$\begin{aligned} \frac{\epsilon _1}{\epsilon _2}=-\nu _{yx}=-\frac{4\mu +8k_1(1-3\kappa )^2 c^2 s^2}{8\mu +8k_1(1-3\kappa )^2c^4} \end{aligned}$$

(52)

And the modulus $E_y$:

$$\begin{aligned} E_y& = \frac{\sigma _{22}}{\epsilon _{22}}=-\frac{(4\mu +8k_1 (1-3\kappa )^2c^2s^2)^2}{8\mu +8k_1(1-3\kappa )^2c^4} \nonumber \\&+\quad (8\mu +8k_1(1-3\kappa )^2s^4)\nonumber \\ E_y& = \frac{6\mu ^2+8k_1\mu (1-3\kappa )^2(1-3c^2s^2)}{\mu +k_1(1-3\kappa )^2c^4} \end{aligned}$$

(53)

The Poisson’s ratio $\nu _{yz}$ is:

$$\begin{aligned}&\epsilon _1+\epsilon _2=-\epsilon _3 \implies 1+\frac{\epsilon _1}{\epsilon _2}=-\frac{\epsilon _3}{\epsilon _2}\nonumber \\&\implies \nu _{yz}=-\frac{\epsilon _3}{\epsilon _2}=\frac{4\mu +8k_1(1-3\kappa )^2 c^2(c^2-s^2)}{8\mu +8k_1(1-3\kappa )^2c^4} \end{aligned}$$

(54)

1.1.3 Pulling in the z direction to determine $E_z, \nu _{zx}, \nu _{zy}$

When the block is subjected to tension in z direction and is free to expand in x, y directions, the stress state becomes:

$$\begin{aligned} \sigma _{11} = 0, \sigma _{22} = 0, \sigma _{33} \ne 0 \end{aligned}$$

(55)

Using $\lambda _2=\frac{1}{\lambda _1\lambda _3}$, Eq. (33) becomes:

$$\begin{aligned} \sigma _{11}& = -p+2\mu \lambda _1^2+4\kappa k_1 E_1 \hbox {exp}[k_2E_1^2] \lambda _1^2 \nonumber \\&+\quad 4(1-3\kappa )k_1E_1\hbox {exp}[k_2E_1^2]c^2\lambda _1^2=0\nonumber \\ \sigma _{22}& = -p+2\mu \frac{1}{\lambda _1^2\lambda _3^2}+4\kappa k_1 E_1 \hbox {exp}[k_2E_1^2] \frac{1}{\lambda _1^2\lambda _3^2} \nonumber \\&\quad +\,4(1-3\kappa ) k_1E_1\hbox {exp}[k_2E_1^2]\frac{s^2}{\lambda _1^2\lambda _3^2}=0\nonumber \\ \sigma _{33}& = -p+2\mu \lambda _3^2+4\kappa k_1E_1\hbox {exp}[k_2E_1^2]\lambda _3^2 \end{aligned}$$

(56)

$$\begin{aligned} \sigma _{33}& = \sigma _{33}-\sigma _{22}=2\mu \left[ \lambda _3^2- \frac{1}{\lambda _1^2\lambda _3^2}\right] \nonumber \\&+\quad \,4\kappa k_1E_1\hbox {exp}[k_2E_1^2\left[ \lambda _3^2-\frac{1}{\lambda _1^2\lambda _3^2}\right] \nonumber \\&\quad -\,4k_1E_1(1-3\kappa )\hbox {exp}[k_2E_1^2]\frac{s^2}{\lambda _1^2\lambda _3^2} \nonumber \\ \sigma _{11}& = \sigma _{11}-\sigma _{22}=2\mu \left[ \lambda _1^2-\frac{1}{\lambda _1^2\lambda _3^2}\right] \nonumber \\&\quad +\,4\kappa k_1E_1\hbox {exp}[k_2E_1^2 \left[ \lambda _1^2-\frac{1}{\lambda _1^2\lambda _3^2}\right] \nonumber \\&\quad +\,4k_1E_1(1-3\kappa )\hbox {exp}[k_2E_1^2]\left[ c^2\lambda _1^2- \frac{s^2}{\lambda _1^2\lambda _3^2}\right] \end{aligned}$$

(57)

Using small deformation approximations:

$$\begin{aligned} \lambda _1& = 1+\epsilon _1, \lambda _1^2=1+2\epsilon _1, \frac{1}{\lambda _1^2}=1-2\epsilon _1\nonumber \\ \lambda _3& = 1+\epsilon _3, \lambda _3^2=1+2\epsilon _3, \frac{1}{\lambda _3^2}=1-2\epsilon _3\nonumber \\ \frac{1}{\lambda _1^2\lambda _3^2}& = 1-2\epsilon _1-2\epsilon _3, \lambda _1^2-\frac{1}{\lambda _1^2\lambda _3^2} \nonumber \\& = 4\epsilon _1+2\epsilon _3, \lambda _3^2-\frac{1}{\lambda _1^2\lambda _3^2}=2\epsilon _1+4\epsilon _3 \nonumber \\ \frac{s^2}{\lambda _1^2\lambda _3^2}& = s^2(1-2\epsilon _1-2\epsilon _3)\nonumber \\ c^2\lambda _1^2-\frac{s^2}{\lambda _1^2\lambda _3^2}& = c^2-s^2+ 2\epsilon _1+2s^2\epsilon _3\nonumber \\ E_1& = (1-3\kappa )(c^2\lambda _1^2+\frac{s^2}{\lambda _1^2 \lambda _3^2}-1) \nonumber \\& = (1-3\kappa )\left[ 2(c^2-s^2)\epsilon _1-2s^2\epsilon _3\right] \end{aligned}$$

(58)

Substituting Eq. (58) into Eq. (57):

$$\begin{aligned} \sigma _{33}& = 2\mu (4\epsilon _3+2\epsilon _1)-4k_1(1-3\kappa )^2 s^2 \nonumber \\&\quad \left[ 2(c^2-s^2)\epsilon _1-2s^2\epsilon _3\right] \nonumber \\ \sigma _{11}& = 2\mu (4\epsilon _1+2\epsilon _3)+4k_1(1-3\kappa )^2 (c^2-s^2) \nonumber \\&\quad \left[ 2(c^2-s^2)\epsilon _1-2s^2\epsilon _3\right] =0 \end{aligned}$$

(59)

Therefore, the Poisson’s ratio $\nu _{zx}$ is:

$$\begin{aligned} \frac{\epsilon _1}{\epsilon _3}=-\nu _{zx}=-\frac{4\mu -8k_1 (1-3\kappa )^2(c^2s^2-s^4)}{8\mu +8k_1(1-3\kappa )^2(c^2-s^2)^2} \end{aligned}$$

(60)

The modulus $E_z$ is:

$$\begin{aligned} E_z=\frac{6\mu ^2+8\mu k_1 (1-3\kappa )^2(1-3c^2s^2)}{\mu +k_1(1-3\kappa )^2(c^2-s^2)^2} \end{aligned}$$

(61)

Poisson’s ratio $\nu _{zy}$

$$\begin{aligned} \nu _{zy}=\frac{4\mu +8k_1(1-3\kappa )^2(c^2-s^2)c^2}{8\mu +8k_1(1-3\kappa )^2(c^2-s^2)^2} \end{aligned}$$

(62)

1.2 Apply simple shear stress $\sigma _{xy}$ to determine shear modulus $G_{xy}$

The deformation gradient and all related quantities to compute stress state in this loading case are:

$$\begin{aligned} F& = \begin{bmatrix} 1 &{}\quad \gamma &{}\quad 0 \\ 0 &{}\quad 1 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 1 \end{bmatrix}, \quad B=\begin{bmatrix} 1+\gamma ^2 &{}\quad \gamma &{}\quad 0 \\ \gamma &{}\quad 1 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 1 \end{bmatrix}, \nonumber \\ C& = \begin{bmatrix} 1 &{}\quad \gamma &{}\quad 0 \\ \gamma &{}\quad 1+\gamma ^2 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 1 \end{bmatrix} \end{aligned}$$

(63)

Other quantities in Eq. (50) for computing Cauchy stresses.

$$\begin{aligned} FL_1& = \begin{bmatrix} 1 &{}\quad \gamma &{}\quad 0 \\ 0 &{}\quad 1 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 1 \end{bmatrix} \begin{bmatrix} c\\ s\\ 0 \end{bmatrix} = \begin{bmatrix} c+s\gamma \\ s\\ 0 \end{bmatrix}, \nonumber \\ FL_2& = \begin{bmatrix} 1 &{}\quad \gamma &{}\quad 0 \\ 0 &{}\quad 1 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 1 \end{bmatrix} \begin{bmatrix} c\\ -s\\ 0 \end{bmatrix} = \begin{bmatrix} c-s\gamma \\ -s\\ 0 \end{bmatrix} \end{aligned}$$

(64)

$$\begin{aligned} I_{41}& = \begin{bmatrix} c&\quad s&\quad 0 \end{bmatrix} \begin{bmatrix} 1 &{}\quad \gamma &{}\quad 0 \\ \gamma &{}\quad 1+\gamma ^2 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 1 \end{bmatrix} \begin{bmatrix} c\\ s\\ 0 \end{bmatrix} \nonumber \\& = 1+2cs\gamma \end{aligned}$$

(65)

$$\begin{aligned} I_{42}& = \begin{bmatrix} c&-s&0 \end{bmatrix} \begin{bmatrix} 1 &{}\quad \gamma &{}\quad 0 \\ \gamma &{}\quad 1+\gamma ^2 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 1 \end{bmatrix} \begin{bmatrix} c\\ -s\\ 0 \end{bmatrix} \nonumber \\& = 1-2cs\gamma \end{aligned}$$

(66)

$$\begin{aligned} E_1& = (1-3\kappa )(2cs\gamma ),\nonumber \\ E_2& = (1-3\kappa )(-2cs\gamma ) \end{aligned}$$

(67)

$$\begin{aligned} FL_1 \otimes FL_1& = \begin{bmatrix} c+s\gamma \\ s \\ 0 \end{bmatrix} \begin{bmatrix} c+s\gamma&\quad s&\quad 0 \end{bmatrix} \nonumber \\& = \begin{bmatrix} (c+s\gamma )^2 &{}\quad s(c+s\gamma ) &{}\quad 0\\ s(c+s\gamma ) &{}\quad s^2 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 0 \end{bmatrix} \end{aligned}$$

(68)

$$\begin{aligned} FL_2 \otimes FL_2& = \begin{bmatrix} c-s\gamma \\ - s \\ 0 \end{bmatrix} \begin{bmatrix} c-s\gamma&\quad - s&\quad 0 \end{bmatrix} \nonumber \\& = \begin{bmatrix} (c-s\gamma )^2 &{}\quad -s(c-s\gamma ) &{}\quad 0\\ -s(c-s\gamma ) &{}\quad s^2 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 0 \end{bmatrix} \end{aligned}$$

(69)

Neglecting high order term of $\gamma$, the shear stress $\sigma _{xy}$ in Eq. (33) becomes:

$$\begin{aligned} \sigma _{xy}=\left[ 2\mu +8k_1(1-3\kappa )^2 c^2 s^2\right] \gamma \end{aligned}$$

(70)

Thus, the shear stress $G_{xy}$ becomes:

$$\begin{aligned} G_{xy}=\frac{\sigma _{xy}}{\gamma }=2\mu +8k_1(1-3\kappa )^2 c^2 s^2 \end{aligned}$$

(71)

With similar approach, the shear moduli $G_{xz}$, $G_{yz}$ can also be determined. For fiber plane (x–y), these shear moduli are equal to the value of isotropic case $G_{yz}=G_{xz}=2\mu$.

Figure 17 illustrates how the longitudinal ($E_x$) and transverse ($E_y$) stiffnesses vary with respect to the fiber angle $\theta$ at two values of fiber stiffness $k_1/\mu _{\mathrm{s}}=5$ and $k_1/\mu _{\mathrm{s}}=2$.

Fig. 17

Equivalent longitudinal and transverse moduli $E_x, E_y$ versus fiber angle $\theta$ for the case of fiber dispersion $\kappa =0$ in (x–y) fiber plane