A unified approach to hyperelastic modelling of plant cell walls

Abstract

The structure of the microfibrils network poses the key challenge to understand not only the mechanical behaviours but also the growth mechanisms of the plant cell wall. A number of competing hypotheses with contradictory assumptions, including the recently proposed contact model and the more conventional tethering model, have been proposed to address the challenge. In the present study, a fibre-reinforced hyperelastic model is proposed aiming to cover both the contact model and the tethering model in a unified mathematical framework. A variety of anisotropies, including those arising from the moiré pattern and the contacts, are taken into account. The formulation is presented in terms of the principal stretches for more direct correlation with the biological experimental observations. The details of the stresses and tangent moduli are provided for readily implementing in finite-element modelling. The simple shear problem is studied to provide the quantitative insight into the long-term biological hypothesis that cell wall growth is driven by shear deformations. The moiré pattern under tensile stretch is demonstrated. The local motion of microfibrils is examined by using rod theory. A corresponding potential growth mode is proposed with the analysis using Burgers vector and torsion in material space.

1. Introduction

The plant cell wall is a multilayer composite in which cellulose microfibrils (CMFs) are playing the role as the fibre reinforcement while the matrix is composed of hemicellulose, pectins, and proteins. A number of publications provided the extensive reviews about the recent progresses [1–9] and the top open questions [10] on the multiscale modelling of the cell wall and its growth from different points of view. Noticeably, both the progresses and open questions focus on the understanding of the cell wall structure, the integration mechanisms of the wall components, and the corresponding mechanical regulations of growth at nano- and larger scale.

As the major load-bearing framework, the structure of the microfibrils network poses the key challenge but also serves as the key clue in the research of the cell wall. A number of hypotheses have been proposed to address the challenge. The different hypotheses, however, may contain the very different, even contradictory, assumptions about the CMFs network.

One of the recently proposed hypotheses of interest is that, the CMFs make a load-bearing network via close physical contacts with one another in bundled regions that potentially function as sites of cell wall loosening and creep [1,2]. The contact, which may be also called a ‘hotspot’ [1], is a thin patch of xyloglucan trapped between the adjacent CMFs bundles so as to play the role of the glue to connect the adjacent CMFs bundles together [3]. This is different from the more conventional tethering model in which the mechanical connections between the CMFs were made by extended xyloglucan chains (part of hemicellulose) that adhered noncovalently to CMF surfaces [11,12].

In the tethering model, the cell wall could be modelled as a composite in which the CMFs as fibres are embedded into and connected via the matrix composed of hemicelluloses and pectins (see Figure 1(a)). Hence the growth of the cell wall was supposed to be regulated by the yielding of matrix even though the CMFs provided the key mechanical stiffness [12–14].

Figure 1.

Two CMFs network models: (a) tethering model in which CMFs are imbedded into and connected by the matrix composing of hemicellulose, pectins and proteins; (b) contact model where CMFs show waviness and are connected directly by some contacts, which is proposed to be geometrically characterised by two assumed families of fibres with the waviness angle $2\mathrm{\delta \theta }$ between them.

By contrast, in the contact model, the direct contacts (connections) between the CMFs may play the key roles in the regulation of the expansive growth of the cell wall as the model suggests that the expansion of the cell wall is mainly constrained by the contacts rather than the matrix surrounding the CMFs [1]. The CMFs may be considered as a interconnected nano-fibre network on its own which, with a weak coupling with the matrix, may independently regulating both the mechanical stiffness and growth of the cell wall (see Figure 1(b)). This is different from the postulate in the tethering model that the CMFs are relative inert and playing a passive role in growth regulation [12].

From the point of view of fibre-reinforced material modelling, in the tethering model, the CMFs at a representative material point of one lamina (layer) of the cell wall may be modelled as one family of the unidirectional fibres imbedded in the matrix. By contrast, in the contact model, the direct connections between the CMFs need to be taken into account. Therefore, one of the questions of interest is if there is a unified mathematical modelling framework to cover both the tethering model and the contact model so that both models may be formulated and validated in a comparable way.

In this study, a finite strain fibre-reinforced hyperelastic model is proposed to address this question. On one hand, the proposed method may cover the tethering model straightforwardly as a material with one-family of reinforced fibres. On the other hand, to reflect the hypothesis of the contact model, rather than one-family of the actual CMFs, it is proposed that there are two assumed families of fibres in the lamina (see Figure 1(b)). The angle, $2\mathrm{\delta \theta }$ , between the two assumed families of fibres represents the waviness of the configurations of the actual CMFs while the directions of the two assumed families of fibres are aligned to the two major directions of the actual CMFs.

It is noted that, by comparison to the tethered model, the contact model may contain additional anisotropies besides the anisotropy due to the presence of physical fibres (CMFs) which is in the scope of the well-established theory of fibre-reinforced materials. This study postulates two classes of additional anisotropies. The first is the anisotropy of moiré pattern which is created by two neighbouring (overlapping) cell wall laminae. Moiré pattern usually is created by a small misalignment of two layers of 2D material, e.g., bilayer graphene [15–18]. In the contact model of cell wall, given that two neighbouring laminae have a finite relative rotation, the small waviness angle $2\mathrm{\delta \theta }$ may play the role of misalignment so as to create moiré pattern, which yields its own anisotropy. Second, the contacts themselves may contain their own anisotropies arising from (1) the shape of the contact as a thin patch of xyloglucan (geometry anisotropy), (2) the arrangement of xyloglucan in the patch (material anisotropy), and (3) the spatial arrangement of a group of contacts (spatial anisotropy).

Taking all the above anisotropies into account, both the tethering model and the contact model may be formulated in a unified mathematical framework. The challenge is to formulate in a way that may be directly correlated with the biological experiment observations, especially the roles of the contacts in the regulation of mechanical and growth behaviours.

Geometrically, a cell wall is a thin-walled structure. Therefore, the stretches of the wall surfaces are the prominent indications of the mechanical behaviour of the cell wall [19], especially in the in vivo experimental measurements. There are a number of hyperelastic models of cell wall reported in literature (e.g. Huang et al. [20,21]) which were mainly adopted and developed from the widely used soft-tissue model proposed by Gasser et al. [22] formulated using the invariants of the Cauchy–Green tensor. For the more direct correlation with the biological experiment observations of the stretches of the cell wall, this study formulates the constitutive model based on the methodology in the principal stretches proposed by Hill [23], Ogden [24], Storåkers [25], and Shariff [26,27]. Noting that the principal stretch formulation of hyperelasticity may have singularity issue when two or three of the principal stretches are identical. This poses the computational difficulty to the numerical modelling like finite-element method (FEM). To address this difficulty, the formulation is presented using the framework proposed by Simo and Taylor [28] so as to provide a method not only for analytical study but also readily to be implemented to the numerical study using FEM. It is worth mentioning that Shariff and co-workers [26,27] have developed the spectral constitutive modelling of transversely isotropic materials which showed concise structure, especially when dealing with analytical problems, and could be considered as the alternative modelling framework.

The outline of the research is summarised here. The fibre-reinforced hyperelastic model taking into the aforementioned variety of anisotropies is proposed in the principal stretches. For the purpose of numerical modelling (e.g. FEM), the formulations of stresses and tangential moduli are discussed in detail. With the unified hyperelastic model formulated, two specific cases, pure shear and tensile stretch, are studied for further understanding the relation between cell wall material and growth. The simple shear problem is analysed, aiming to give the quantitative insight into one of the key biological hypotheses that cell wall growth is driven/achieved by shear deformations [1,13,14]. On the other hand, tensile stretch serves the case to examine the following specific novel topics: (1) the moiré pattern under tensile stretch, (2) CMFs’ local motion when a shared contact is loosened and the adjacent CMFs are separated by tensile stretch, which may be modelled by using the rod theory at microscopic scale; and consequently (3) the potential growth mode if the previously connected CMFs can separate locally, for which Burgers vector and torsion in material space are introduced for the analysis of the incompatible growth configuration.

2. Kinematics of the model

The kinematics may be described in the framework of the conventional continuum theory of a closed system as the present study is interested in the elastic behaviour of the cell wall. The mapping $\chi :{\mathcal{B}}_0\to {\mathcal{B}}_t$ is defined as a deformation which maps the reference configuration, ${\mathcal{B}}_0$ , onto the current configuration, ${\mathcal{B}}_t$ , so that the current position of a material point X is written as $\mathit{x}=\chi (\mathit{X},t)\in {\mathcal{B}}_t$ , where $t_0$ is dropped from the expression for conciseness. The deformation gradient is defined as $\mathit{F}(\mathit{X},t)=\frac{\mathrm{\partial \chi }}{\partial \mathit{X}}=\frac{\partial {\chi }_i}{\partial {\mathit{X}}_I}{\mathit{e}}_i\otimes {\mathit{E}}_I$ with its Jacobian denoted as $J=\det \mathit{F}>0$ . Whereby ${\mathit{E}}_I$ ( $I=1,2,3$ ) and ${\mathit{e}}_i$ ( $i=1,2,3$ ) are the fixed orthonormal bases in ${\mathcal{B}}_0$ and ${\mathcal{B}}_t$ , respectively.

2.1. Kinematics of CMFs network

Similar to engineering composite material, the cell wall has a multilayered laminae structure. We shall first discuss the kinematics of CMFs in a single lamina which is applied to all individual layers. Second, the additional periodic structures arising from the moiré pattern created by two overlapping layers are discussed. Third, we discuss the description of the contacts taking into account its spatial density and anisotropy.

2.1.1. Kinematics of CMFs in the individual lamina of cell wall

For the tethering model, the CMFs in Figure 1(a) may be geometrically represented by a family of fibres assigned a specific reference direction at the material point X on ${\mathcal{B}}_0$ . On the other hand, for the contact model shown in Figure 1(b), the CMFs may be considered as two families of fibres with their directions different by the waviness angle, $2\mathrm{\delta \theta }$ . In this way, both tethering model and contact model may be modelled in a unified kinematical framework.

Let a unit vector, ${\mathit{a}}_{I0}(\mathit{X})$ , denote the reference fibre direction of the $I^{\mathrm{th}}$ -family $(I\in \{1,\dots ,n_f\})$ of fibres attached to the solid particle at the material point X on ${\mathcal{B}}_0$ . We are particularly interested in the case of $n_f=2$ for contact model as well as $n_f=1$ for tethering model. But let us assume that $n_f$ can be any finite integral number for the general purpose of modelling. By using the perfectly bonding assumption, the vector ${\mathit{a}}_{I0}(\mathit{X})$ is mapped to ${\mathcal{B}}_t$ by a push-forward as

$${\mathit{a}}_I(\mathit{X},t)=\mathit{F}(\mathit{X},t){\mathit{a}}_{I0}(\mathit{X}).$$

(1)

Correspondingly, we may define the structural tensors, ${\mathit{A}}_{I0}={\mathit{a}}_{I0}\otimes {\mathit{a}}_{I0}$ and ${\mathit{a}}_{I0}\otimes {\mathit{a}}_{J0}(I\ne J)$ on ${\mathcal{B}}_0$ and ${\mathit{A}}_I={\mathit{a}}_I\otimes {\mathit{a}}_I$ and ${\mathit{a}}_I\otimes {\mathit{a}}_J(I\ne J)$ on ${\mathcal{B}}_0$ , respectively.

2.1.2. Interaction of CMFs of two laminae of cell wall

We consider the bilayer periodic structure by using the moiré pattern theory which has been widely applied in the study of 2D materials like bilayer graphene and quantum materials [15–18]. This study suggests that the waviness angle $2\mathrm{\delta \theta }$ may play the role of misalignment to create the additional periodicity of moiré pattern when two layers of cell wall laminae are taken into account.

Let us consider two typical overlapping laminae in the cell wall (see Figure 2). It is assumed that both layers have the identical structure of the contact model as shown in Figure 1(b) with the waviness angle $2\mathrm{\delta \theta }$ . The second layer is rotated rigidly by a finite angle $\mathrm{\Lambda }(\gg \mathrm{\delta \theta })$ with respect to the first layer. Hence, given the averaged fibre direction, ${\mathit{a}}_0$ , in the first layer, the set of the two-families $(n_f=2)$ of CMFs is

$${\mathcal{B}}_0^{\{1\}}=\left\{{\mathit{a}}_{10}^{\{1\}}=\hat{\mathit{R}}(-\mathrm{\delta \theta }){\mathit{a}}_0,{\mathit{a}}_{20}^{\{1\}}=\hat{\mathit{R}}(\mathrm{\delta \theta }){\mathit{a}}_0\right\},$$

(2)

Figure 2.

The forming of moiré pattern by two overlapped cell wall laminae: the two layers have a finite relative rotation by an angle Λ while the half waviness angle $\mathrm{\delta \theta }(\ll \mathrm{\Lambda })$ in the contact model plays the role of the misalignment of two layers of 2D material to create the moiré pattern.

while in the second layer, CMFs are represented as a set,

$${\mathcal{B}}_0^{\{2\}}=\left\{{\mathit{a}}_{10}^{\{2\}}=\hat{\mathit{R}}(\mathrm{\Lambda }){\mathit{a}}_{10}^{\{1\}}=\hat{\mathit{R}}(\mathrm{\Lambda }-\mathrm{\delta \theta }){\mathit{a}}_0,{\mathit{a}}_{20}^{\{2\}}=\hat{\mathit{R}}(\mathrm{\Lambda }){\mathit{a}}_{20}^{\{1\}}=\hat{\mathit{R}}(\mathrm{\Lambda }+\mathrm{\delta \theta }){\mathit{a}}_0\right\}.$$

(3)

Whereby $\hat{\mathit{R}}(\ast )$ is an in-plane 2D rotation by an angle ∗, the superscripts $\{1\}$ and $\{2\}$ refer to the first and second layers, respectively. We may simply write the above two sets as

$${\mathcal{B}}_0^{\{J\}}=\left\{{\mathit{a}}_{I0}^{\{J\}}|I\in \{1,2\}\right\}\mathrm{\forall }J=1,2.$$

It is clear that ${\mathcal{B}}_0^{\{2\}}$ can be obtained from ${\mathcal{B}}_0^{\{1\}}$ by imposing the finite rotation $\hat{\mathit{R}}(\mathrm{\Lambda })$ . This does not meet the requirement of small angle misalignment to create moiré pattern. Noting that the waviness angle $2\mathrm{\delta \theta }$ is small by comparison to the interlayer rotation angle Λ. Hence, instead of the individual layer structure of CMFs characterised by ${\mathcal{B}}_0^{\{J\}}$ ( $J=1,2$ ), we may consider such two sets,

$${\mathcal{B}}_0^{<I>=\left\{{\mathit{a}}_{I0}^{\{J\}}|J\in \{1,2\}\right\}\mathrm{\forall }I=1,2,}$$

for the formulation of the moiré pattern. ${\mathcal{B}}_0^{<2>}$ is obtained from ${\mathcal{B}}_0^{<1>}$ by imposing a small rotation $\hat{\mathit{R}}(2\mathrm{\delta \theta })$ , which geometrically leads to the moiré pattern. The corresponding moiré vectors, which characterise the moiré pattern, may be defined by a double reciprocal operation as [18]:

$${\mathit{v}}_{0K}=\Re \left(\Re \left({\mathit{a}}_{10}^{\{K\}}\right)-\Re \left({\mathit{a}}_{20}^{\{K\}}\right)\right)\mathrm{\forall }K=1,2,$$

(4)

where $\Re (\ast )$ stands for the reciprocal vector of ∗. Note that the above definition is different from the one proposed by Escudero et al. [17] where only one reciprocal operation, i.e., the difference of the reciprocal vectors, was calculated.

The explicit expression of moiré vectors ${\mathit{v}}_{0K}$ in equation (4) are

$${\mathit{v}}_{0K}={(\mathbf{1}-\hat{\mathit{R}}(2\mathrm{\delta \theta }))}^{-\mathrm{T}}{\mathit{a}}_{10}^{\{K\}}\mathrm{\forall }K=1,2,$$

(5)

where the components of the unit tensor 1 has the value of Kronecker delta, i.e., ${[\mathbf{1}]}_{\mathrm{ij}}={\delta }_{\mathrm{ij}}$ . Then the normalised moiré vectors ${\overline{\mathit{v}}}_{0K}$ are defined as

$${\overline{\mathit{v}}}_{0K}=\frac{{\mathit{v}}_{0K}}{||{\mathit{v}}_{0K}||}\mathrm{\forall }K=1,2,$$

(6)

where the symbol $||\ast ||\stackrel{def}{=}\sqrt{\ast \cdot \ast }$ stands for the Euclidean norm of a vector ∗. Explicitly,

$$||{\mathit{v}}_{0K}||=\sqrt{{\mathit{a}}_{10}^{\{K\}}\cdot {\mathit{H}}_0{\mathit{a}}_{10}^{\{K\}}},$$

(7)

in which

$${\mathit{H}}_0={(2\mathbf{1}-(\hat{\mathit{R}}(2\mathrm{\delta \theta })+{\hat{\mathit{R}}}^{\mathrm{T}}(2\mathrm{\delta \theta })))}^{-1}.$$

(8)

The example of moiré pattern and moiré vectors are presented in Figure 3.

Figure 3.

The moiré pattern formed by two overlapped cell wall laminae where the contact model is used with the half waviness angle $\mathrm{\delta \theta }=5^{\mathrm{\circ }}$ . The two layers have a $\mathrm{\Lambda }=60^{\mathrm{\circ }}$ relative rotation. The white arrows indicate the moiré vectors.

The sets ${\mathcal{B}}_0^{<I>}$ in the reference configuration ${\mathcal{B}}_0$ may be pushed forward to the current configuration ${\mathcal{B}}_t$ as

$${\mathcal{B}}_t^{<I>=\left\{{\overline{\mathit{a}}}_I^{\{J\}}|J\in \{1,2\}\right\}\mathrm{\forall }I=1,2,}$$

where ${\overline{\mathit{a}}}_I^{\{J\}}={\mathit{a}}_I^{\{J\}}∕||{\mathit{a}}_I^{\{J\}}||$ while ${\mathit{a}}_I^{\{J\}}=F{\mathit{a}}_{I0}^{\{J\}}$ ( $\mathrm{\forall }I,J=1,2$ in which I for fibres and J for layers). Here it is assumed that the deformation gradients F in both layers are the same. It is possible to be generalised to the situation where different layers are strained differently (see e.g. Escudero et al. [17]).

Then the moiré vectors on ${\mathcal{B}}_t$ are:

$${\mathit{v}}_K=\Re \left(\Re \left({\overline{\mathit{a}}}_1^{\{K\}}\right)-\Re \left({\overline{\mathit{a}}}_2^{\{K\}}\right)\right)\mathrm{\forall }K=1,2.$$

(9)

The explicit expression of moiré vectors ${\mathit{v}}_K$ are

$${\mathit{v}}_K=\mathit{F}{(\mathbf{1}-\hat{\mathit{R}}(2\mathrm{\delta \theta }))}^{-\mathrm{T}}{\mathit{a}}_{10}^{\{K\}}∕||{\mathit{a}}_1^{\{K\}}||\mathrm{\forall }K=1,2,$$

(10)

where $||{\mathit{a}}_1^{\{K\}}||=||\mathit{F}{\mathit{a}}_{10}^{\{K\}}||$ . The normalised moiré vectors ${\overline{\mathit{v}}}_K$ are defined as

$${\overline{\mathit{v}}}_K=\frac{{\mathit{v}}_K}{||{\mathit{v}}_K||}\mathrm{\forall }K=1,2,$$

(11)

where

$$||{\mathit{v}}_K||=\frac{1}{||{\mathit{a}}_1^{\{K\}}||}\sqrt{{\mathit{a}}_{10}^{\{K\}}\cdot \mathit{H}{\mathit{a}}_{10}^{\{K\}}},$$

(12)

in which

$$\mathit{H}={(\mathbf{1}-\hat{\mathit{R}}(2\mathrm{\delta \theta }))}^{-1}\mathit{C}{(\mathbf{1}-\hat{\mathit{R}}(2\mathrm{\delta \theta }))}^{-\mathrm{T}}.$$

(13)

We have the following two relations between the moiré vectors in ${\mathcal{B}}_t$ and in ${\mathcal{B}}_0$ for $\mathrm{\forall }K=1,2$ . The first is

$${\mathit{v}}_K=\frac{\mathit{F}{\mathit{v}}_{0K}}{||{\mathit{a}}_1^{\{K\}}||}=\frac{\mathit{F}{\mathit{v}}_{0K}}{||\mathit{F}{\mathit{a}}_{10}^{\{K\}}||},$$

(14)

which is in the same form of the relation between the normalised fibre direction in ${\mathcal{B}}_t$ and the unit fibre in ${\mathcal{B}}_0$ , i.e., ${\overline{\mathit{a}}}_1^{\{K\}}={\mathit{a}}_1^{\{K\}}∕||{\mathit{a}}_1^{\{K\}}||=\mathit{F}{\mathit{a}}_{01}^{\{K\}}∕||\mathit{F}{\mathit{a}}_{10}^{\{K\}}||$ . The second one is

$${\overline{\mathit{v}}}_K=\frac{\mathit{F}{\mathit{v}}_{0K}}{||{\mathit{a}}_1^{\{K\}}||||{\mathit{v}}_K||}=\frac{||{\mathit{v}}_{0K}||}{||{\mathit{a}}_1^{\{K\}}||||{\mathit{v}}_K||}\mathit{F}{\overline{\mathit{v}}}_{0K}.$$

(15)

Analogues to ${\mathit{A}}_{I0}$ and ${\mathit{A}}_I$ , the structural tensors, ${\mathit{V}}_{0K}={\overline{\mathit{v}}}_{0K}\otimes {\overline{\mathit{v}}}_{0K}$ and ${\mathit{V}}_K={\mathit{v}}_K\otimes {\mathit{v}}_K$ ( $\mathrm{\forall }K=1,2$ ), may be defined to represent the effects of moiré pattern, which reflects the interlayer effects.

The moiré angles, ${\omega }_0$ in ${\mathcal{B}}_0$ and ω in ${\mathcal{B}}_t$ , can be calculated as

$${\omega }_0=\arccos ({\overline{\mathit{v}}}_{01}\cdot {\overline{\mathit{v}}}_{02})\mathrm{and}\omega =\arccos ({\overline{\mathit{v}}}_1\cdot {\overline{\mathit{v}}}_2),$$

(16)

respectively. It can be proven that

$${\omega }_0=\arccos ({\mathit{a}}_{10}^{\{K\}}\cdot {\mathit{a}}_{20}^{\{K\}}).$$

(17)

2.1.3. Description of the contacts

In the reference configuration, the contacts, which physically are the thin patches of xyloglucan between CMFs bundles, are considered as a discrete point set, ${\mathcal{P}}_0\subset {\mathcal{B}}_0$ , embedded into the matrix while connecting the CMFs into a network. Each point in the set ${\mathcal{P}}_0$ is defined by a representative material point of a contact, e.g., the mass centre of the corresponding xyloglucan patch. We assumed that the point set is perfectly bonded to the matrix during elastic deformations so that the mapping χ defines the current set ${\mathcal{P}}_t=\chi ({\mathcal{P}}_0)\subset {\mathcal{B}}_t$ .

It is proposed that the point set of contacts ${\mathcal{P}}_0$ and ${\mathcal{P}}_t$ may be characterised by a scalar, ζ, and a structural tensor, κ . The scalar ζ may be assigned as (1) the concentration of the contact points (number per unit volume) in a representative elementary volume (REV) in ${\mathcal{B}}_0$ , (2) the characteristic length between two contact points (200nm reported by Zhang et al. [4]), and (3) the averaged ratio of the contact size to the length between two contact points. Noting that the contacts can be loosened during growth which may lead to the decrease of number of contacts or contact size and the increase of the length between two contact points. Hence, ζ may evolve during growth accordingly.

On the other hand, the structural tensor κ may represent the anisotropic characteristics of the contact as a patch of xyloglucan taking into account: (1) the geometric anisotropy of the shape of the patch, (2) the material anisotropy of the arrangement of the xyloglucan in the patch, and (3) the spatial anisotropy of the arrangement of the contacts in REV which may be referred to the mathematical method developed for analysing the crystal structures [29,30].

2.1.4. Summary of the anisotropies of CMFs network

The anisotropies of CMFs network introduced in the above discussion are summarised in Table 1.

Table 1.

Anisotropies of CMFs network.

Type	Layer(s)	Anisotropies	Structural tensors
1	Monolayer	Fibre direction(s)	${\mathit{A}}_{I0}$ and ${\mathit{A}}_I$
2	Bilayer	Moiré vectors	${\mathit{V}}_{0K}$ and ${\mathit{V}}_K$
3	Monolayer	Contacts:	κ
	/Multilayers	(1) Shape of patch,
		(2) Material arrangement in a patch,
		(3) Spatial arrangement of patches.

Mathematically, the structural tensors ${\mathit{V}}_{0K}$ and ${\mathit{V}}_K$ of moiré pattern are similar to those of fibres, ${\mathit{A}}_{I0}$ and ${\mathit{A}}_I$ . The formulation using in the latter may be similarly applied to the former. Hence we omit the details of the formulation related to ${\mathit{V}}_{0K}$ and ${\mathit{V}}_K$ in the following discussion where they can be obtained by simply transforming the counterpart of the formulation of ${\mathit{A}}_{I0}$ and ${\mathit{A}}_I$ .

On the other hand, the details of formulation related to κ are presented as they are not in the identical forms as the formulation of ${\mathit{A}}_{I0}$ and ${\mathit{A}}_I$ .

2.2. Kinematics in principal stretches

The polar decompositions of F are

$$\mathit{F}=\mathit{RU}=\mathit{VR},$$

(18)

where R is a rotation tensor while U and $\mathit{V}$ are the right and left stretch tensors, respectively. Accordingly, the right and left Cauchy-Green tensors, C and b , are defined as

$$\mathit{C}={\mathit{F}}^{\mathrm{T}}\mathit{F}\equiv {\mathit{U}}^2\mathrm{and}\mathit{b}=\mathit{F}{\mathit{F}}^{\mathrm{T}}\equiv {\mathit{V}}^2,$$

(19)

respectively.

By using the spectral theorem for symmetric positive definite tensors, C may be expressed in terms of the principal stretches, ${\lambda }_i$ (i=1,2,3), as [31]

$$\mathit{C}=\sum _{i=1}^3{\lambda }_i^2{\mathit{N}}_i\otimes {\mathit{N}}_i,$$

(20)

where $||{\mathit{N}}_i||=1$ for the Euclidean norm ||•|| of the vector ${\mathit{N}}_i$ , and ${\lambda }_i^2>0$ ( $i=1,2,3$ ) are the roots of the characteristic polynomial

$$p({\lambda }_i^2)\stackrel{def}{=}-{\lambda }_i^6+I_1{\lambda }_i^4-I_2{\lambda }_i^2+I_3=0,$$

(21)

in which $I_i$ ( $i=1,2,3$ ) are the three principal invariants of C ,

$$I_1=tr\mathit{C},I_2=\frac{1}{2}[I_1^2-tr{\mathit{C}}^2],I_3=\det \mathit{C}=J^2,$$

(22)

respectively.

Accordingly, equation (19) indicates that

$$\mathit{U}=\sum _{i=1}^3{\lambda }_i{\mathit{N}}_i\otimes {\mathit{N}}_i.$$

(23)

And the spectral decompositions of F and R may be expressed as

$$\mathit{F}=\sum _{i=1}^3{\lambda }_i{\mathit{n}}_i\otimes {\mathit{N}}_i\mathrm{and}\mathit{R}=\sum _{i=1}^3{\mathit{n}}_i\otimes {\mathit{N}}_i,$$

(24)

respectively, where ${\mathit{n}}_i$ is defined by a transformation,

$${\lambda }_i{\mathit{n}}_i=\mathit{F}{\mathit{N}}_i.$$

(25)

Consequently, from equation (19), b is expressed in terms of ${\mathit{n}}_i$ as

$$\mathit{b}=\sum _{i=1}^3{\lambda }_i^2{\mathit{n}}_i\otimes {\mathit{n}}_i.$$

(26)

For convenience, a structural tensor, ${\mathit{M}}_i$ , may be defined as

$${\mathit{M}}_i={\lambda }_i^{-2}{\mathit{N}}_i\otimes {\mathit{N}}_i\equiv {\mathit{F}}^{-1}({\mathit{n}}_i\otimes {\mathit{n}}_i){\mathit{F}}^{-\mathrm{T}}.$$

(27)

Then, for example, C can be expressed as $\mathit{C}={\mathrm{\sum }}_{i=1}^3{\lambda }_i^4{\mathit{M}}_i$ .

Noting that, although ${\lambda }_i^2$ need to be solved from equation (21), ${\mathit{N}}_i\otimes {\mathit{N}}_i$ and ${\mathit{n}}_i\otimes {\mathit{n}}_i$ are not necessary to be solved directly by the spectral decomposition of C and b . Instead, there are the closed form solutions [28,32],

$${\mathit{N}}_i\otimes {\mathit{N}}_i=\frac{{\lambda }_i^2}{D_i}\left[\mathit{C}-(I_1-{\lambda }_i^2)\mathbf{1}+I_3{\lambda }_i^{-2}{\mathit{C}}^{-1}\right]$$

(28)

and

$${\mathit{n}}_i\otimes {\mathit{n}}_i=\frac{1}{D_i}\left[{\mathit{b}}^2-(I_1-{\lambda }_i^2)\mathit{b}+I_3{\lambda }_i^{-2}\mathbf{1}\right],$$

(29)

if and only if $D_i\ne 0$ , where

$$D_i=2{\lambda }_i^4-I_1{\lambda }_i^2+I_3{\lambda }_i^{-2}=({\lambda }_i^2-{\lambda }_j^2)({\lambda }_i^2-{\lambda }_k^2)\mathrm{\forall }i\ne j\ne k.$$

(30)

This indicates the singularity issue, $D_i=0$ , when either two of the principal stretches are equal (e.g. ${\lambda }_1={\lambda }_2\ne {\lambda }_3$ ) or all of them are identical (i.e. ${\lambda }_1={\lambda }_2={\lambda }_3$ ). According to this observation, we shall address all three typical cases, Case 1: ${\lambda }_1\ne {\lambda }_2\ne {\lambda }_3$ , Case 2: ${\lambda }_1={\lambda }_2\ne {\lambda }_3$ , and Case 3: ${\lambda }_1={\lambda }_2={\lambda }_3$ , separately, for the constitutive modelling [28].

2.3. Generalised strain measures

In either $\mathrm{in}\mathrm{vivo}$ or $\mathrm{in}\mathrm{vitro}$ experimental measurements, the principal stretches ${\lambda }_i(i=1,2,3)$ of the cell wall surface are among the variables in principle observable and measurable so as to be effectively used in the plant science community [19]. This, as aforementioned, is one of the observations motivates this study to formulate the elasticity using the principal stretches.

In order to effectively correlate the potential experimental data with, say, a power law of any principal stretch, the generalised strain measures are introduced as [31],

$${\mathit{E}}^{(m)}=\{\begin{array}{cc}\frac{1}{m}[{\mathit{C}}^{m∕2}-\mathbf{1}]\hfill & \mathrm{if}m\ne 0\hfill \\ \frac{1}{2}\ln \mathit{C}\hfill & \mathrm{if}m=0\hfill \end{array},$$

(31)

where, by using equation (20), the expression of ${\mathit{C}}^{m∕2}$ for Case 1-Case 3 of the different states of ${\lambda }_i$ (i=1,2,3) are

$${\mathit{C}}^{m∕2}=\{\begin{array}{cc}\sum _{i=1}^3{\lambda }_i^m{\mathit{N}}_i\otimes {\mathit{N}}_i\hfill & \mathrm{for}\mathrm{Case}1:{\lambda }_1\ne {\lambda }_2\ne {\lambda }_3\hfill \\ {\lambda }_1^m\mathbf{1}+({\lambda }_3^m-{\lambda }_1^m){\mathit{N}}_3\otimes {\mathit{N}}_3\hfill & \mathrm{for}\mathrm{Case}2:{\lambda }_1={\lambda }_2\ne {\lambda }_3\hfill \\ {\lambda }_1^m\mathbf{1}\hfill & \mathrm{for}\mathrm{Case}3:{\lambda }_1={\lambda }_2={\lambda }_3\hfill \end{array}.$$

(32)

2.4. The derivative relations between ${\lambda }_i$ and $\mathit{C}$

Besides the tensorial relation (20), the derivative relations between ${\lambda }_i(i=1,2,3)$ and C (and its invariants $I_i$ ( $i=1,2,3$ )) play the key role in the hyperelastic formulation (and numerical study) using the principal stretches ${\lambda }_i(i=1,2,3)$ [31]. The relations are provided here for the completeness.

By assuming ${\lambda }_1\ne {\lambda }_2\ne {\lambda }_3$ , equation (21) indicates that ${\lambda }_i$ ( $\mathrm{\forall }i=1,2,3$ ) is a function of the invariants of C , say

$${\lambda }_i={\hat{\lambda }}_i(I_1,I_2,I_3)={\overline{\lambda }}_i(\mathit{C})\mathrm{\forall }i=1,2,3.$$

(33)

To avoid over complexity in the expressions, we simply use ${\lambda }_i$ in the following discussion when referring to the functions ${\hat{\lambda }}_i(I_1,I_2,I_3)$ and ${\overline{\lambda }}_i(\mathit{C})$ .

First, the derivatives of ${\lambda }_i$ with respect to $I_j(\mathrm{\forall }j=1,2,3)$ are [31]

$$\frac{\partial {\lambda }_i}{\partial I_1}=\frac{1}{2}\frac{{\lambda }_i^3}{D_i},\frac{\partial {\lambda }_i}{\partial I_2}=-\frac{1}{2}\frac{{\lambda }_i}{D_i},\frac{\partial {\lambda }_i}{\partial I_3}=\frac{1}{2}\frac{{\lambda }_i^{-1}}{D_i},$$

(34)

which requires that $D_i\ne 0$ .

Second, the calculation of $\mathrm{d}{\lambda }_i∕\mathrm{d}\mathit{C}$ is given as follows.

2.4.1. Case 1: ${\lambda }_1\ne {\lambda }_2\ne {\lambda }_3$

Differentiation of equation (20) in both sides reads

$$\mathrm{d}\mathit{C}=\sum _{i=1}^3(2{\lambda }_i\mathrm{d}{\lambda }_i{\mathit{N}}_i\otimes {\mathit{N}}_i+{\lambda }_i^2\mathrm{d}{\mathit{N}}_i\otimes {\mathit{N}}_i+{\lambda }_i^2{\mathit{N}}_i\otimes \mathrm{d}{\mathit{N}}_i).$$

(35)

Since $\mathrm{d}{\mathit{N}}_i\cdot {\mathit{N}}_i=0$ due to $||{\mathit{N}}_i||=1$ , the double contraction of equation (35) by ${\mathit{N}}_i\otimes {\mathit{N}}_i$ yields

$$\frac{\mathrm{d}{\lambda }_i}{\mathrm{d}\mathit{C}}=\frac{1}{2{\lambda }_i}{\mathit{N}}_i\otimes {\mathit{N}}_i=\frac{{\lambda }_i}{2}{\mathit{M}}_i\mathrm{\forall }i=1,2,3.$$

(36)

2.4.2. Case 2: ${\lambda }_1={\lambda }_2\ne {\lambda }_3$

In this situation,

$$\mathit{C}={\lambda }_1^2\left(\sum _{i=1}^2{\mathit{N}}_i\otimes {\mathit{N}}_i\right)+{\lambda }_3^2{\mathit{N}}_3\otimes {\mathit{N}}_3.$$

(37)

Hence, using equation (35) we have

$$\frac{\mathrm{d}{\lambda }_1}{\mathrm{d}\mathit{C}}=\frac{1}{2{\lambda }_1}\left(\sum _{i=1}^2{\mathit{N}}_i\otimes {\mathit{N}}_i\right)=\frac{{\lambda }_1}{2}({\mathit{C}}^{-1}-{\mathit{M}}_3)$$

(38)

and

$$\frac{\mathrm{d}{\lambda }_3}{\mathrm{d}\mathit{C}}=\frac{1}{2{\lambda }_3}({\mathit{N}}_3\otimes {\mathit{N}}_3)=\frac{{\lambda }_3}{2}{\mathit{M}}_3,$$

(39)

where the relation,

$${\mathit{C}}^{-1}-{\mathit{M}}_3={\lambda }_1^{-2}\sum _{i=1}^2{\mathit{N}}_i\otimes {\mathit{N}}_i,$$

is derived from equations (25) and (27).

2.4.3. Case 3: ${\lambda }_1={\lambda }_2={\lambda }_3$

We write equation (20) as

$$\mathit{C}={\lambda }_1^2\mathbf{1}.$$

(40)

Hence $\mathrm{d}\mathit{C}=2{\lambda }_1\mathrm{d}{\lambda }_1\mathbf{1}$ , which means

$$\frac{\mathrm{d}{\lambda }_{\mathbf{1}}}{\mathrm{d}\mathit{C}}=\frac{1}{2{\lambda }_1}\mathbf{1}=\frac{{\lambda }_1}{2}{\mathit{C}}^{-1}.$$

(41)

3. Constitutive theory

3.1. Free energy functions of the cell wall as a fibre-reinforced material

Let us consider the free energy functions (i.e. free energy densities per unit volume) on ${\mathcal{B}}_0$ as

$$W(\mathit{X},\mathit{F},{\mathit{A}}_{I0},{\mathit{V}}_{0K},\zeta ,\mathit{\kappa })=\overline{W}(\mathit{X},\mathit{C},{\mathit{A}}_{I0},{\mathit{V}}_{0K},\zeta ,\mathit{\kappa })$$

(42)

$$=Ŵ(\mathit{X},{\lambda }_1,{\lambda }_2,{\lambda }_3)+{\overline{W}}_f(\mathit{X},\mathit{C},{\mathit{A}}_{I0})+{\overline{W}}_m(\mathit{X},\mathit{C},{\mathit{V}}_{K0})+W_c(\mathit{X},\mathit{C},\zeta ,\mathit{\kappa }),$$

in which Ŵ, ${\overline{W}}_f$ , ${\overline{W}}_m$ and $W_c$ represent the matrix, fibre, moiré vector, and contact contributions to W (or $\overline{W}$ ), respectively, $I\in \{1,\dots ,n_f\}$ refers to all families of fibres in the REV of the material point X , and $K\in \{1,\dots ,n_l-1\}$ represents the $K^{\mathrm{th}}$ -pair of neighbouring (overlapping) layers given $n_l$ the total number of cell wall lamimea.

For conciseness, it is useful to introduce a function ${\beta }_i$ ,

$${\beta }_i=\frac{\mathrm{\partial Ŵ}}{\partial {\lambda }_i}{\lambda }_i(\mathrm{no}\mathrm{summation}),$$

(43)

and calculate its derivative with respect to the stretches as

$$\frac{\partial {\beta }_i}{\partial {\lambda }_j}=\frac{{\partial }^2Ŵ}{\partial {\lambda }_i\partial {\lambda }_j}{\lambda }_i+\frac{\mathrm{\partial Ŵ}}{\partial {\lambda }_i}{\delta }_{\mathrm{ij}}.$$

(44)

3.1.1. Specification of the function $Ŵ$

In order to model the matrix composed of hemicellulose and pectins, we consider two features of the cell wall matrix. First, the recent reports referred the pectins as a hydrogel gel [1,5]. Second, the plant cell wall subjected to dehydration may still offer the useful insight into its microstructure [33]. The two features suggest that a matrix model capable of modelling a wide range of volumetric deformation from nearly incompressible to highly compressible may provide a broad support for exploiting and interpreting experimental data.

Therefore, a specification of the free energy function Ŵ of the wall matrix in equation (42) is adopted from the isotropic compressible model proposed by Storåkers [25]:

$$\begin{array}{c}Ŵ(\mathit{X},{\lambda }_1,{\lambda }_2,{\lambda }_3)\equiv {Ŵ}_0(\mathit{X},{\mathit{C}}^{m_k∕2})\hfill \\ =\sum _{k=1}^{n_m}\frac{c_k}{m_k}\left(\left(I_1({\mathit{C}}^{m_k∕2})-3\right)+\frac{1}{n}(I_3^{-n}({\mathit{C}}^{m_k∕2})-1)\right),\hfill \end{array}$$

(45)

where $c_k$ , $m_k$ ( $\mathrm{\forall }k\in \{1,\dots ,n_m\}$ for a certain integer $n_m$ ), and n are the material parameters. The first and third invariants of ${\mathit{C}}^{m_k∕2}$ (see equation (32) for $m=m_k$ ) can be expressed in terms of ${\lambda }_i$ as

$$I_1({\mathit{C}}^{m_k∕2})=\sum _{i=1}^3{\lambda }_i^{m_k},I_3({\mathit{C}}^{m_k∕2})={({\lambda }_1{\lambda }_2{\lambda }_3)}^{m_k}=J^{m_k}.$$

Hence, Ŵ in equation (45) may be expressed explicitly in terms of ${\lambda }_i$ as

$$Ŵ={Ŵ}_0={Ŵ}_1(\mathit{X},{\lambda }_1,{\lambda }_2,{\lambda }_3)=\sum _{k=1}^{n_m}\frac{c_k}{m_k}\left(\sum _{i=1}^3{\lambda }_i^{m_k}-3+\frac{1}{n}(J^{-m_{k}n}-1)\right).$$

(46)

Accordingly, the specifications of ${\beta }_i$ and $\partial {\beta }_i∕\partial {\lambda }_j$ in equations (43) and (44) are

$${\beta }_i=\frac{\mathrm{\partial Ŵ}}{\partial {\lambda }_i}{\lambda }_i=\frac{\partial {Ŵ}_1}{\partial {\lambda }_i}{\lambda }_i=\sum _{k=1}^{n_m}{c}_k\left({\lambda }_i^{m_k}-J^{-m_{k}n}\right),$$

(47)

and

$$\frac{\partial {\beta }_i}{\partial {\lambda }_j}=\sum _{k=1}^{n_m}{c}_k{m}_k\left({\lambda }_i^{m_k-1}{\delta }_{\mathrm{ij}}+n{J}^{-m_{k}n}{\lambda }_j^{-1}\right).$$

3.1.2. Specification of the functions ${\overline{W}}_f$ and ${\overline{W}}_m$

The free energy function ${\overline{W}}_f(\mathit{X},\mathit{C},{\mathit{A}}_{I0})$ of the fibres is supposed in a form as

$$\begin{array}{c}{\overline{W}}_f(\mathit{X},\mathit{C},{\mathit{A}}_{I0})={Ŵ}_f(\mathit{X},I_{I04}^{(n_k)},I_{\mathrm{IJ}04}^{(n_k)})\hfill \\ =\frac{1}{2}\sum _{I=1}^{n_f}\sum _{k=1}^{n_n}{\overline{c}}_I^{(k)}{\left(I_{I04}^{(n_k)}-1\right)}^2+\frac{1}{2}\sum _{I,J=1(I<J)}^{n_f}\sum _{k=1}^{n_n}{\tilde{c}}_{\mathrm{IJ}}^{(k)}{\left(I_{\mathrm{IJ}04}^{(n_k)}-I_{\mathrm{IJ}04}^{(n_k)}{|}_{\mathit{C}=\mathbf{1}}\right)}^2\hfill \\ (\mathrm{\forall }I,J\in \{1,\dots ,n_f\}),\hfill \end{array}$$

(48)

where ${\overline{c}}_I^{(k)}$ , ${\tilde{c}}_{\mathrm{IJ}}^{(k)}$ , and $n_k$ ( $k\in \{1,\dots ,n_n\}$ for a certain integer $n_n$ ) are the material parameters, $I_{I04}^{(n_k)}$ is the invariant of the $I^{\mathrm{th}}$ -family of fibres while $I_{\mathrm{IJ}04}^{(n_k)}$ the coupling invariant of the $I^{\mathrm{th}}$ -family and $J^{\mathrm{th}}$ -family of fibres. In consistency with the discussion of the matrix deformation in terms of the principle stretches, the expressions of $I_{I04}^{(n_k)}$ and $I_{\mathrm{IJ}04}^{(n_k)}$ are discussed for Case 1-Case 3 of the different states of stretches, respectively.

Case 1: ${\lambda }_1\ne {\lambda }_2\ne {\lambda }_3$

Using equations (20) and (27), it is defined that

$$I_{I04}^{(n_k)}\stackrel{def}{=}tr({\mathit{C}}^{n_k∕2}{\mathit{A}}_{I0})=\sum _{i=1}^3{\lambda }_i^{n_k+2}tr({\mathit{M}}_i{\mathit{A}}_{I0}),$$

(49)

and

$$I_{\mathrm{IJ}04}^{(n_k)}\stackrel{def}{=}tr({\mathit{C}}^{n_k∕2}[{\mathit{a}}_{I0}\otimes {\mathit{a}}_{J0}])=\sum _{i=1}^3{\lambda }_i^{n_k+2}tr({\mathit{M}}_i[{\mathit{a}}_{I0}\otimes {\mathit{a}}_{J0}]),$$

(50)

respectively. Noting that the two specifications of $I_{I04}^{(n_k)}$ are equivalent to the usual invariants, $I_{I04}$ and $I_{I06}$ , of C ,

$$\begin{array}{c}{I}_{I04}^{(2)}=I_{I04}=tr(\mathit{C}{\mathit{A}}_{I0})\mathrm{if}{n}_k=2,\hfill \\ I_{I04}^{(4)}=I_{I06}=tr({\mathit{C}}^2{\mathit{A}}_{I0})\mathrm{if}{n}_k=4.\hfill \end{array}$$

Case 2: ${\lambda }_1={\lambda }_2\ne {\lambda }_3$

In this case, $I_{I04}^{(n_k)}$ and $I_{\mathrm{IJ}04}^{(n_k)}$ are expressed as

$$I_{I04}^{(n_k)}={\lambda }_1^{n_k+2}tr(({\mathit{C}}^{-1}-{\mathit{M}}_3){\mathit{A}}_{I0})+{\lambda }_3^{n_k+2}tr({\mathit{M}}_3{\mathit{A}}_{I0})$$

(51)

and

$$I_{\mathrm{IJ}04}^{(n_k)}={\lambda }_1^{n_k+2}tr(({\mathit{C}}^{-1}-{\mathit{M}}_3)[{\mathit{a}}_{I0}\otimes {\mathit{a}}_{J0}])+{\lambda }_3^{n_k+2}tr({\mathit{M}}_3[{\mathit{a}}_{I0}\otimes {\mathit{a}}_{J0}]),$$

(52)

respectively.

Case 3: ${\lambda }_1={\lambda }_2={\lambda }_3$

$I_{I04}^{(n_k)}$ and $I_{\mathrm{IJ}04}^{(n_k)}$ in this case have simple forms:

$$I_{I04}^{(n_k)}={\lambda }_1^{n_k}tr{\mathit{A}}_{I0}$$

(53)

and

$$I_{\mathrm{IJ}04}^{(n_k)}={\lambda }_1^{n_k}tr[{\mathit{a}}_{I0}\otimes {\mathit{a}}_{J0}].$$

(54)

The function ${\overline{W}}_m(\mathit{X},\mathit{C},{\mathit{V}}_{K0})$ can be defined in the similar form as ${\overline{W}}_f(\mathit{X},\mathit{C},{\mathit{A}}_{I0})$ with ${\mathit{A}}_{I0}$ replaced by ${\mathit{V}}_{K0}$ . The details are omitted here as aforementioned.

3.1.3. Specification of the function $W_c$

It is supposed that ${\overline{W}}_c$ has such a general form

$$W_c(\mathit{X},\mathit{C},\zeta ,\mathit{\kappa })=\hat{\varphi }(\zeta ){Ŵ}_c(\mathit{X},{\lambda }_1,{\lambda }_2,{\lambda }_3)+\overline{\varphi }(\zeta ){\overline{W}}_c(\mathit{X},\mathit{C},\mathit{\kappa }),$$

(55)

where $\hat{\varphi }(\zeta )$ and $\overline{\varphi }(\zeta )$ are two functions which represent the effect of wall loosening. It is expected that $\hat{\varphi }(\zeta )$ and $\overline{\varphi }(\zeta )$ decrease when wall is loosened via the loosening of the contacts, which reduces $W_c$ accordingly.

The two functions, ${Ŵ}_c(\mathit{X},{\lambda }_1,{\lambda }_2,{\lambda }_3)$ and ${\overline{W}}_c(\mathit{X},\mathit{C},\mathit{\kappa })$ , may have the similar form as $Ŵ(\mathit{X},{\lambda }_1,{\lambda }_2,{\lambda }_3)$ and ${\overline{W}}_f(\mathit{X},\mathit{C},{\mathit{A}}_{I0})$ , respectively. Therefore, it is assumed that

$${Ŵ}_c(\mathit{X},{\lambda }_1,{\lambda }_2,{\lambda }_3)=\sum _{k=1}^{{\hat{n}}_m}\frac{{ĉ}_k}{{\hat{m}}_k}\left(\left(I_1({\mathit{C}}^{{\hat{m}}_k∕2})-3\right)+\frac{1}{\hat{n}}(I_3^{-\hat{n}}({\mathit{C}}^{{\hat{m}}_k∕2})-1)\right),$$

(56)

where ${ĉ}_k$ , ${\hat{m}}_k$ ( $\mathrm{\forall }k\in \{1,\dots ,{\hat{n}}_m\}$ for a certain integer ${\hat{n}}_m$ ), and $\hat{n}$ are the material parameters, and

$${\overline{W}}_c(\mathit{X},\mathit{C},\mathit{\kappa })={\hat{\overline{W}}}_c(\mathit{X},I_{\mathit{\kappa }}^{({\hat{n}}_k)})=\frac{1}{2}\sum _{k=1}^{{\hat{n}}_n}{\overline{c}}_{\mathit{\kappa }}^{(k)}{\left(I_{\mathit{\kappa }}^{({\hat{n}}_k)}-1\right)}^2,$$

(57)

where $I_{\mathit{\kappa }}^{({\hat{n}}_k)}$ ( $\mathrm{\forall }k\in \{1,\dots ,{\hat{n}}_n\}$ for a certain integer ${\hat{n}}_n$ ) are defined as:

Case 1: ${\lambda }_1\ne {\lambda }_2\ne {\lambda }_3$

$$I_{\mathit{\kappa }}^{({\hat{n}}_k)}\stackrel{def}{=}tr({\mathit{C}}^{{\hat{n}}_k∕2}\mathit{\kappa })=\sum _{i=1}^3{\lambda }_i^{{\hat{n}}_k+2}tr({\mathit{M}}_i\mathit{\kappa }),$$

(58)

Case 2: ${\lambda }_1={\lambda }_2\ne {\lambda }_3$

$$I_{\mathit{\kappa }}^{({\hat{n}}_k)}={\lambda }_1^{{\hat{n}}_k+2}tr(({\mathit{C}}^{-1}-{\mathit{M}}_3)\mathit{\kappa })+{\lambda }_3^{{\hat{n}}_k+2}tr({\mathit{M}}_3\mathit{\kappa }),$$

(59)

Case 3: ${\lambda }_1={\lambda }_2={\lambda }_3$

$$I_{\mathit{\kappa }}^{({\hat{n}}_k)}={\lambda }_1^{{\hat{n}}_k}tr\mathit{\kappa }.$$

(60)

3.2. Stresses

The second Piola–Kirchhoff (PK2) stress, S , on ${\mathcal{B}}_0$ and Cauchy stress, σ , on ${\mathcal{B}}_t$ are obtained from the hyperelastic relations,

$$\mathit{S}=2\frac{\partial \overline{W}}{\partial \mathit{C}}\mathrm{and}\mathit{\sigma }=J^{-1}\mathit{F}\mathit{S}{\mathit{F}}^{\mathrm{T}}=J^{-1}\mathit{F}\frac{\mathrm{\partial W}}{\partial \mathit{F}},$$

(61)

respectively. According to the decomposition of the total free energy function in equation (42), the total stresses S and σ may be additively decomposed as

$$\mathit{S}=\hat{\mathit{S}}+{\hat{\mathit{S}}}_f+{\hat{\mathit{S}}}_m+{\hat{\mathit{S}}}_c\mathrm{and}\mathit{\sigma }=\hat{\mathit{\sigma }}+{\hat{\mathit{\sigma }}}_f+{\hat{\mathit{\sigma }}}_m+{\hat{\mathit{\sigma }}}_c,$$

(62)

where

$$\hat{\mathit{S}}=2\frac{\mathrm{\partial Ŵ}}{\partial \mathit{C}},{\hat{\mathit{S}}}_f=2\frac{\partial {\overline{W}}_f}{\partial \mathit{C}},{\hat{\mathit{S}}}_m=2\frac{\partial {\overline{W}}_m}{\partial \mathit{C}},\mathrm{and}{\hat{\mathit{S}}}_c=2\frac{\partial W_c}{\partial \mathit{C}}$$

(63)

are the matrix part, fibre part, moiré part, and contact part, respectively, of S , while $\hat{\mathit{\sigma }}$ , ${\hat{\mathit{\sigma }}}_f$ , ${\hat{\mathit{\sigma }}}_m$ , and ${\hat{\mathit{\sigma }}}_c$ are the counterparts corresponding to σ .

3.2.1. Matrix part of stresses

Using $\partial {\lambda }_i∕\partial \mathit{C}$ in equations (36), (38), (39), and (41), equation (63)₁ indicates that

$$\begin{array}{c}\hat{\mathit{S}}=2\frac{\mathrm{\partial Ŵ}(\mathit{X},{\lambda }_1,{\lambda }_2,{\lambda }_3)}{\partial {\lambda }_i}\frac{\partial {\lambda }_i}{\partial \mathit{C}}\hfill \\ =\{\begin{array}{cc}\sum _{i=1}^3{\beta }_i{\mathit{M}}_i\hfill & \mathrm{for}\mathrm{Case}1:{\lambda }_1\ne {\lambda }_2\ne {\lambda }_3\hfill \\ {\beta }_1{\mathit{C}}^{-1}+({\beta }_3-{\beta }_1){\mathit{M}}_3\hfill & \mathrm{for}\mathrm{Case}2:{\lambda }_1={\lambda }_2\ne {\lambda }_3\hfill \\ {\beta }_1{\mathit{C}}^{-1}\hfill & \mathrm{for}\mathrm{Case}3:{\lambda }_1={\lambda }_2={\lambda }_3\hfill \end{array}.\hfill \end{array}$$

(64)

Based on equation (27), we have the relations $\mathit{F}{\mathit{M}}_i{\mathit{F}}^{\mathrm{T}}={\mathit{n}}_i\otimes {\mathit{n}}_i$ . Hence, Cauchy stress $\hat{\mathit{\sigma }}$ is calculated from equation (64) by using the same transform as in equation (61)₂,

$$\begin{array}{c}\hat{\mathit{\sigma }}=J^{-1}\mathit{F}\hat{\mathit{S}}{\mathit{F}}^{\mathrm{T}}\hfill \\ =J^{-1}\{\begin{array}{cc}\hfill \sum _{i=1}^3{\beta }_i{\mathit{n}}_i\otimes {\mathit{n}}_i\hfill & \hfill \mathrm{for}\mathrm{Case}1:{\lambda }_1\ne {\lambda }_2\ne {\lambda }_3\hfill \\ \hfill {\beta }_1\mathbf{1}+\left({\beta }_3-{\beta }_1\right){\mathit{n}}_i\otimes {\mathit{n}}_i\hfill & \hfill \mathrm{for}\mathrm{Case}2:{\lambda }_1={\lambda }_2\ne {\lambda }_3\hfill \\ \hfill {\beta }_1\mathbf{1}\hfill & \hfill \mathrm{for}\mathrm{Case}3:{\lambda }_1={\lambda }_2={\lambda }_3\hfill \end{array}.\hfill \end{array}$$

(65)

3.2.2. Fibre and moiré parts of stresses

Using equation (63)₂ together with equation (48), we write the fibre part of the PK2 stress as

$$\begin{array}{c}{\hat{\mathit{S}}}_f=2\frac{\partial {Ŵ}_f(\mathit{X},I_{I04}^{(n_k)},I_{\mathrm{IJ}04}^{(n_k)})}{\partial \mathit{C}}=2\sum _{I=1}^{n_f}\sum _{k=1}^{n_n}\frac{\partial {Ŵ}_f(\mathit{X},I_{I04}^{(n_k)},I_{\mathrm{IJ}04}^{(n_k)})}{\partial I_{I04}^{(n_k)}}\frac{\partial I_{I04}^{(n_k)}}{\partial \mathit{C}}\hfill \\ +2\sum _{I,J=1(I<J)}^{n_f}\sum _{k=1}^{n_n}\frac{\partial {Ŵ}_f(\mathit{X},I_{I04}^{(n_k)},I_{\mathrm{IJ}04}^{(n_k)})}{\partial I_{\mathrm{IJ}04}^{(n_k)}}\frac{\partial I_{\mathrm{IJ}04}^{(n_k)}}{\partial \mathit{C}}.\hfill \end{array}$$

(66)

To be consistent with the formulation of the matrix part of stresses in equation (64), ${\hat{\mathit{S}}}_f$ may be formulated for the different states of stretches as described in Case 1–Case 3. For this purpose, we shall discuss the expressions of the two terms, $\partial I_{I04}^{(n_k)}∕\partial \mathit{C}$ and $\partial I_{\mathrm{IJ}04}^{(n_k)}∕\partial \mathit{C}$ , in the right-hand-side (r.h.s) of equation (66) for Case 1-Case 3.

Case 1: ${\lambda }_1\ne {\lambda }_2\ne {\lambda }_3$

First, the term $\partial I_{I04}^{(n_k)}∕\partial \mathit{C}$ in the r.h.s of equation (66) is calculated by using equation (58) as

$$\frac{\partial I_{I04}^{(n_k)}}{\partial \mathit{C}}=\sum _{i=1}^{3}tr({\mathit{M}}_i{\mathit{A}}_{I0})(n_k+2){\lambda }_i^{n_k+1}\frac{\partial {\lambda }_i}{\partial \mathit{C}}+\sum _{i=1}^3{\lambda }_i^{n_k+2}{\mathit{A}}_{I0}:\frac{\partial {\mathit{M}}_i}{\partial \mathit{C}},$$

(67)

where $\partial {\lambda }_i∕\partial \mathit{C}$ is provided in equation (36), the calculation of the term $\partial {\mathit{M}}_i∕\partial \mathit{C}$ is provided in Appendix 1 for completeness (see equation (161) for the expression).

Second, the coupling term $\partial I_{\mathrm{IJ}04}^{(n_k)}∕\partial \mathit{C}$ is calculated with equation (50),

$$\frac{\partial I_{\mathrm{IJ}04}^{(n_k)}}{\partial \mathit{C}}=\sum _{i=1}^{3}tr({\mathit{M}}_i[{\mathit{a}}_{I0}\otimes {\mathit{a}}_{J0}])(n_k+2){\lambda }_i^{n_k+1}\frac{\partial {\lambda }_i}{\partial \mathit{C}}+\sum _{i=1}^3{\lambda }_i^{n_k+2}[{\mathit{a}}_{I0}\otimes {\mathit{a}}_{J0}]:\frac{\partial {\mathit{M}}_i}{\partial \mathit{C}}.$$

(68)

Case 2: ${\lambda }_1={\lambda }_2\ne {\lambda }_3$

Similarly, the term $\partial I_{I04}^{(n_k)}∕\partial \mathit{C}$ in this case is calculated by using equation (59) as

$$\begin{array}{c}\frac{\partial I_{I04}^{(n_k)}}{\partial \mathit{C}}=(n_k+2){\lambda }_1^{n_k+1}tr(({\mathit{C}}^{-1}-{\mathit{M}}_3){\mathit{A}}_{I0})\frac{\partial {\lambda }_1}{\partial \mathit{C}}+(n_k+2){\lambda }_3^{n_k+1}tr({\mathit{M}}_3{\mathit{A}}_{I0})\frac{\partial {\lambda }_3}{\partial \mathit{C}}\hfill \\ -{\lambda }_1^{n_k+2}{\mathit{A}}_{I0}:{\mathbb{I}}_{{\mathit{C}}^{-1}}+({\lambda }_3^{n_k+2}-{\lambda }_1^{n_k+2}){\mathit{A}}_{I0}:\frac{\partial {\mathit{M}}_3}{\partial \mathit{C}},\hfill \end{array}$$

(69)

where we already use the relation,

$$\frac{\partial {\mathit{C}}^{-1}}{\partial \mathit{C}}=-{\mathbb{I}}_{{\mathit{C}}^{-1}},\mathrm{in}\mathrm{which}{[{\mathbb{I}}_{{\mathit{C}}^{-1}}]}_{\mathrm{ijkl}}=\frac{1}{2}(C_{\mathrm{ik}}^{-1}{C}_{\mathrm{jl}}^{-1}+C_{\mathrm{il}}^{-1}{C}_{\mathrm{jk}}^{-1}),$$

(70)

and $\partial {\lambda }_1∕\partial \mathit{C}$ and $\partial {\lambda }_3∕\partial \mathit{C}$ are provided in equations (38) and (39), respectively.

For the coupling term $\partial I_{\mathrm{IJ}04}^{(n_k)}∕\partial \mathit{C}$ , equation (52) is used to calculate that

$$\begin{array}{c}\frac{\partial I_{\mathrm{IJ}04}^{(n_k)}}{\partial \mathit{C}}=(n_k+2){\lambda }_1^{n_k+1}tr(({\mathit{C}}^{-1}-{\mathit{M}}_3)[{\mathit{a}}_{I0}\otimes {\mathit{a}}_{J0}])\frac{\partial {\lambda }_1}{\partial \mathit{C}}\hfill \\ +(n_k+2){\lambda }_3^{n_k+1}tr({\mathit{M}}_3[{\mathit{a}}_{I0}\otimes {\mathit{a}}_{J0}])\frac{\partial {\lambda }_3}{\partial \mathit{C}}\hfill \\ -{\lambda }_1^{n_k+2}[{\mathit{a}}_{I0}\otimes {\mathit{a}}_{J0}]:{\mathbb{I}}_{{\mathit{C}}^{-1}}+({\lambda }_3^{n_k+2}-{\lambda }_1^{n_k+2})[{\mathit{a}}_{I0}\otimes {\mathit{a}}_{J0}]:\frac{\partial {\mathit{M}}_3}{\partial \mathit{C}}.\hfill \end{array}$$

(71)

Case 3: ${\lambda }_1={\lambda }_2={\lambda }_3$

In this case, both terms, $\partial I_{I04}^{(n_k)}∕\partial \mathit{C}$ and $\partial I_{\mathrm{IJ}04}^{(n_k)}∕\partial \mathit{C}$ , are calculated straightforwardly by using equations (54) and (60), respectively, as

$$\frac{\partial I_{I04}^{(n_k)}}{\partial \mathit{C}}=n_k{\lambda }_1^{n_k-1}tr{\mathit{A}}_{I0}\frac{\partial {\lambda }_1}{\partial \mathit{C}}$$

(72)

and

$$\frac{\partial I_{\mathrm{IJ}04}^{(n_k)}}{\partial \mathit{C}}=n_k{\lambda }_1^{n_k-1}tr[{\mathit{a}}_{I0}\otimes {\mathit{a}}_{J0}]\frac{\partial {\lambda }_1}{\partial \mathit{C}},$$

(73)

where $\partial {\lambda }_1∕\partial \mathit{C}$ can be referred to equation (41).

${\hat{\mathit{S}}}_m$ can be calculated in the similar way as ${\hat{\mathit{S}}}_f$ if the free energy function ${\overline{W}}_m(\mathit{X},\mathit{C},{\mathit{V}}_{K0})$ is defined in the similar form as ${\overline{W}}_f(\mathit{X},\mathit{C},{\mathit{A}}_{I0})$ .

3.2.3. Contact part of stresses

Using equations (55) and (57), equation (63)₄ indicates that

$${\hat{\mathit{S}}}_c=2\hat{\varphi }(\zeta )\frac{\partial {Ŵ}_c(\mathit{X},{\lambda }_1,{\lambda }_2,{\lambda }_3)}{\partial {\lambda }_i}\frac{\partial {\lambda }_i}{\partial \mathit{C}}+2\overline{\varphi }(\zeta )\frac{\partial {\hat{\overline{W}}}_c(\mathit{X},I_{\mathit{\kappa }}^{({\hat{n}}_k)})}{\partial \mathit{C}}.$$

(74)

The first term in r.h.s of equation (74) is calculated in a form similar to equation (64):

$$2\frac{\partial {Ŵ}_c(\mathit{X},{\lambda }_1,{\lambda }_2,{\lambda }_3)}{\partial {\lambda }_i}\frac{\partial {\lambda }_i}{\partial \mathit{C}}=\{\begin{array}{cc}\sum _{i=1}^3{\hat{\beta }}_i{\mathit{M}}_i\hfill & \mathrm{for}\mathrm{Case}1:{\lambda }_1\ne {\lambda }_2\ne {\lambda }_3\hfill \\ {\hat{\beta }}_1{\mathit{C}}^{-1}+({\hat{\beta }}_3-{\hat{\beta }}_1){\mathit{M}}_3\hfill & \mathrm{for}\mathrm{Case}2:{\lambda }_1={\lambda }_2\ne {\lambda }_3\hfill \\ {\hat{\beta }}_1{\mathit{C}}^{-1}\hfill & \mathrm{for}\mathrm{Case}3:{\lambda }_1={\lambda }_2={\lambda }_3\hfill \end{array}.$$

(75)

in which, similar to equation (47),

$${\hat{\beta }}_i=\frac{\partial {Ŵ}_c}{\partial {\lambda }_i}{\lambda }_i=\sum _{k=1}^{{\hat{n}}_m}{ĉ}_k\left({\lambda }_i^{{\hat{m}}_k}-J^{-{\hat{m}}_k\hat{n}}\right).$$

(76)

The second term in r.h.s of equation (74) is calculated in a form similar to equation (66) as

$$2\frac{\partial {\hat{\overline{W}}}_c(\mathit{X},I_{\mathit{\kappa }}^{({\hat{n}}_k)})}{\partial \mathit{C}}=2\sum _{k=1}^{n_n}\frac{\partial {\hat{\overline{W}}}_c(\mathit{X},I_{\mathit{\kappa }}^{({\hat{n}}_k)})}{\partial I_{\mathit{\kappa }}^{({\hat{n}}_k)}}\frac{\partial I_{\mathit{\kappa }}^{({\hat{n}}_k)}}{\partial \mathit{C}}$$

(77)

The term $\frac{\partial I_{\mathit{\kappa }}^{({\hat{n}}_k)}}{\partial \mathit{C}}$ is calculated for three cases:

Case 1: ${\lambda }_1\ne {\lambda }_2\ne {\lambda }_3$

Similar to equation (67), we have

$$\frac{\partial I_{\mathit{\kappa }}^{({\hat{n}}_k)}}{\partial \mathit{C}}=\sum _{i=1}^{3}tr({\mathit{M}}_i\mathit{\kappa })({\hat{n}}_k+2){\lambda }_i^{{\hat{n}}_k+1}\frac{\partial {\lambda }_i}{\partial \mathit{C}}+\sum _{i=1}^3{\lambda }_i^{{\hat{n}}_k+2}\mathit{\kappa }:\frac{\partial {\mathit{M}}_i}{\partial \mathit{C}},$$

(78)

Case 2: ${\lambda }_1={\lambda }_2\ne {\lambda }_3$

In this case, similar to equation (69),

$$\begin{array}{c}\frac{\partial I_{\mathit{\kappa }}^{({\hat{n}}_k)}}{\partial \mathit{C}}=({\hat{n}}_k+2){\lambda }_1^{{\hat{n}}_k+1}tr(({\mathit{C}}^{-1}-{\mathit{M}}_3)\mathit{\kappa })\frac{\partial {\lambda }_1}{\partial \mathit{C}}+({\hat{n}}_k+2){\lambda }_3^{{\hat{n}}_k+1}tr({\mathit{M}}_3\mathit{\kappa })\frac{\partial {\lambda }_3}{\partial \mathit{C}}\hfill \\ -{\lambda }_1^{{\hat{n}}_k+2}\mathit{\kappa }:{\mathbb{I}}_{{\mathit{C}}^{-1}}+({\lambda }_3^{{\hat{n}}_k+2}-{\lambda }_1^{{\hat{n}}_k+2})\mathit{\kappa }:\frac{\partial {\mathit{M}}_3}{\partial \mathit{C}},\hfill \end{array}$$

(79)

Case 3: ${\lambda }_1={\lambda }_2={\lambda }_3$

Similar to equation (72), the expression is

$$\frac{\partial I_{\mathit{\kappa }}^{({\hat{n}}_k)}}{\partial \mathit{C}}={\hat{n}}_k{\lambda }_1^{{\hat{n}}_k-1}tr\mathit{\kappa }\frac{\partial {\lambda }_1}{\partial \mathit{C}}.$$

(80)

Acorrdingly, Cauchy stress ${\hat{\mathit{\sigma }}}_c$ in equation (62) is calculated from equation (74) by using the transform, ${\hat{\mathit{\sigma }}}_c=J^{-1}\mathit{F}{\hat{\mathit{S}}}_c{\mathit{F}}^{\mathrm{T}}$ .

3.3. Tangent moduli

Besides stresses, the tangent moduli play the key role in numerical modelling, e.g., FEM. The tangent moduli on ${\mathcal{B}}_0$ are

$$\mathrm{ℂ}=4\frac{{\partial }^2\overline{W}}{\partial \mathit{C}\partial \mathit{C}},$$

(81)

while the tangent moduli c on ${\mathcal{B}}_t$ are obtained by using the push-forward of ℂ from ${\mathcal{B}}_0$ to ${\mathcal{B}}_t$ as

$${[\mathbf{c}]}_{\mathrm{ijkl}}=F_{\mathrm{iI}}{F}_{\mathrm{jJ}}{F}_{\mathrm{kK}}{F}_{\mathrm{lL}}{[\mathrm{ℂ}]}_{\mathrm{IJKL}}.$$

(82)

According to the decomposition of free energy density (42) or the decomposition of stresses (62), the total tangent moduli ℂ and c may be additively decomposed as

$$\mathrm{ℂ}=\widehat{\mathrm{ℂ}}+{\widehat{\mathrm{ℂ}}}_f+{\widehat{\mathrm{ℂ}}}_m+{\widehat{\mathrm{ℂ}}}_c\mathrm{and}\mathbf{c}=\hat{\mathbf{c}}+{\hat{\mathbf{c}}}_f+{\hat{\mathbf{c}}}_m+{\hat{\mathbf{c}}}_c,$$

(83)

where

$$\begin{array}{c}\widehat{\mathrm{ℂ}}=4\frac{{\partial }^2Ŵ}{\partial \mathit{C}\partial \mathit{C}}=2\frac{\partial \hat{\mathit{S}}}{\partial \mathit{C}},{\widehat{\mathrm{ℂ}}}_f=2\frac{{\partial }^2{\overline{W}}_f}{\partial \mathit{C}\partial \mathit{C}}=2\frac{\partial {\hat{\mathit{S}}}_f}{\partial \mathit{C}},{\widehat{\mathrm{ℂ}}}_m=2\frac{{\partial }^2{\overline{W}}_m}{\partial \mathit{C}\partial \mathit{C}}=2\frac{\partial {\hat{\mathit{S}}}_m}{\partial \mathit{C}},\hfill \\ {\widehat{\mathrm{ℂ}}}_c=2\frac{{\partial }^2{W}_c}{\partial \mathit{C}\partial \mathit{C}}=2\frac{\partial {\hat{\mathit{S}}}_c}{\partial \mathit{C}}.\hfill \end{array}$$

(84)

are the matrix part, fibre part, moiré part, and contact part, respectively, of ℂ, while $\hat{\mathbf{c}}$ , ${\hat{\mathbf{c}}}_f$ , ${\hat{\mathbf{c}}}_m$ , and ${\hat{\mathbf{c}}}_c$ are the counterparts corresponding to c.

The detailed calculation of tangent moduli ℂ and c is provided in Appendix 2. The results can be implemented in FEM software, e.g., Abaqus^®, by transferring into user’s material subroutine (UMAT) properly.

4. Simple shear of cell wall lamina

One of the prevailing hypotheses of cell wall growth mechanisms is that shear deformations plays the key role, biologically and physically, to achieve the expansive growth of the cell wall [1].

Based on this hypothesis, Boudaoud [13] and Dumais et al. [14] proposed the growth models exactly analogues to the conventional isotropic $J_2$ -plasticity theory [34]. Huang et al. [20,21] developed the more sophisticated finite strain growth models analogues to the anisotropic viscoplasticity [35]. Typically, cell walls were modelled as the thin-walled cylindrical or dome structures by those researchers for the full scale analysis. Especially, with the implementation of the mechanical and growth models in FEM, it is feasible to model various cell wall behaviours including those in the typical mechanical testings like inflation [20,21].

However, analysis of growth at cell wall structural scale mainly provides the solution of wall morphogenesis, which is simply the output of the growth hypothesis. As the shear-induced growth hypothesis was qualitatively presented in biological literature as the collective mechanisms at the material scale of cell wall (e.g. sliding of CMFs) (see Cosgrove [1] for example), it is of great interest to formulate and examine the wall material behaviours under shear deformation for quantitatively representing the biological hypothesis.

Noting that cell wall is a multilayer (laminae) composite structure in which each layer (lamina) can be modelled as a fibre-reinforced material as illustrated in Figure 1. Hence, rather than wall laminae structure, let us consider one representative wall lamina corresponding to Figure 1. Based on the preceding discussion about the hyperelastic modelling of the plant cell wall, we formulate the elastic relation between the simply shear and von Mises stress (equivalent to $J_2$ invariant) in one representative lamina. This closed-form relation may provide the explicitly quantitative insight into the existing hypothesis of shear-induced cell wall growth which, to the author’s best knowledge, has not been reported before in literature. The result serves the purpose to quantitatively represent and validate the shear-induced growth hypothesis without involving wall structural level analysis. The discussion here may be extended to the anisotropic yielding criteria beyond $J_2$ -criterion [20,35] and inelastic analysis taking into account growth law explicitly in future work.

4.1. Kinematics of the simple shear deformation

Let us consider the simple shear deformation [26,27,31,36],

$$\mathit{F}=\left[\begin{array}{ccc}\hfill 1\hfill & \hfill \gamma \hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \end{array}\right],$$

where γ is the in-plane shear of one lamina of the cell wall placed on the ${\mathit{E}}_1\wedge {\mathit{E}}_2$ plane. Then, according to equation (22), the three principal invariants of C are specified as

$$I_1=tr\mathit{C}=3+{\gamma }^2,I_2=\frac{1}{2}[I_1^2-tr{\mathit{C}}^2]=3+{\gamma }^2,I_3=\det \mathit{C}=J^2=1.$$

(85)

The specification of equation (21) reads

$$p({\lambda }_i^2)=-{\lambda }_i^6+(3+{\gamma }^2){\lambda }_i^4-(3+{\gamma }^2){\lambda }_i^2+1=0.$$

(86)

Solving equation (86) gives the principal stretches:

$${\lambda }_{1,2}^2=\frac{(2+{\gamma }^2)\pm \gamma \sqrt{4+{\gamma }^2}}{2},{\lambda }_3^2=1.$$

(87)

By using equation (30), $D_i$ ( $i=1,2,3$ ) are specified as

$$D_1={\gamma }^2(1+{\lambda }_1^2),D_2={\gamma }^2(1+{\lambda }_2^2),D_3=-{\gamma }^2.$$

(88)

The specifications of ${\mathit{M}}_i$ and ${\mathit{M}}_i$ -related entities in the case of simple shear are presented in Appendix 3. Using equation (173) in Appendix 3, we can identify that

$${\mathit{N}}_1=\frac{\gamma }{\sqrt{D_1}}\left\{\begin{array}{c}\hfill 1\hfill \\ \hfill {\lambda }_1\hfill \\ \hfill 0\hfill \end{array}\right\},{\mathit{N}}_2=\frac{\gamma }{\sqrt{D_2}}\left\{\begin{array}{c}\hfill 1\hfill \\ \hfill -{\lambda }_2\hfill \\ \hfill 0\hfill \end{array}\right\},{\mathit{N}}_3=\left\{\begin{array}{c}\hfill 0\hfill \\ \hfill 0\hfill \\ \hfill 1\hfill \end{array}\right\}.$$

(89)

It is straightforward to show that ${\mathit{N}}_1\cdot {\mathit{N}}_2=0.$ Noting that equation (173) is invariant under the change of the sign of principal directions, i.e., ${\mathit{N}}_i\to -{\mathit{N}}_i$ . Moreover, since $\gamma ∕\sqrt{D_1}=\sqrt{1+{\lambda }_1^2}$ and $\gamma ∕\sqrt{D_2}=\sqrt{1+{\lambda }_2^2}$ , the expression (89) does not have singularity issue when γ is approaching zero.

By using the relation (25), i.e., ${\mathit{n}}_i={\lambda }_i^{-1}\mathit{F}{\mathit{N}}_i$ , the spatial vectors, ${\mathit{n}}_i$ $(i=1,2,3)$ , are calculated from equation (89):

$${\mathit{n}}_1=\frac{\gamma }{{\lambda }_1\sqrt{D_1}}\left\{\begin{array}{c}\hfill 1+\gamma {\lambda }_1\hfill \\ \hfill {\lambda }_1\hfill \\ \hfill 0\hfill \end{array}\right\},{\mathit{n}}_2=\frac{\gamma }{{\lambda }_2\sqrt{D_2}}\left\{\begin{array}{c}\hfill 1-\gamma {\lambda }_2\hfill \\ \hfill -{\lambda }_2\hfill \\ \hfill 0\hfill \end{array}\right\},{\mathit{n}}_3=\left\{\begin{array}{c}\hfill 0\hfill \\ \hfill 0\hfill \\ \hfill 1\hfill \end{array}\right\}.$$

(90)

F in equation (24) is invariant under the change of the sign of principal directions, ${\mathit{N}}_i\to -{\mathit{N}}_i$ and ${\mathit{n}}_i\to -{\mathit{n}}_i$ . Equations (87), (89) and (90) fully determine the specification of the spectral decompositions of F . Hence the kinematics in principal stretches is specified.

4.2. Stresses

As of interest is a single lamina, the moiré part and contact part stresses, ${\hat{\mathit{S}}}_m$ and ${\hat{\mathit{S}}}_c$ , are ignored in this section. With the kinematics specified, the contributions of the matrix and fibres to the PK2 stresses, i.e., $\hat{\mathit{S}}$ and ${\hat{\mathit{S}}}_f$ in the decomposition (62), may be calculated from the general formulation, (64) and (66), respectively. Then the total PK2 stress, S , can be obtained as stated in equation (62). The process to calculate the specification of S is provided in Appendix 4.

Consequently, with the help of the relation $\mathit{F}{\mathit{M}}_i{\mathit{F}}^{\mathrm{T}}={\mathit{n}}_i\otimes {\mathit{n}}_i$ (see equation (27)), Cauchy stress σ is calculated from S in equation (179) using the relation (61)₂ as

$$\begin{array}{c}\mathit{\sigma }=J^{-1}\mathit{F}\mathit{S}{\mathit{F}}^{\mathrm{T}}=\sum _{i=1}^3{\beta }_i{J}^{-1}{\mathit{n}}_i\otimes {\mathit{n}}_i+2J^{-1}\sum _{I=1}^{n_f}\sum _{k=1}^{n_n}\{\frac{\partial {Ŵ}_f(\mathit{X},I_{I04}^{(n_k)},I_{\mathrm{IJ}04}^{(n_k)})}{\partial I_{I04}^{(n_k)}}\hfill \\ \left\{\sum _{i=1}^3{\lambda }_i^{n_k+2}\left\{\left(\frac{n_k+2}{2}tr({\mathit{M}}_i{\mathit{A}}_{I0})-\frac{I_3{\lambda }_i^{-2}}{D_i}tr({\mathit{C}}^{-1}{\mathit{A}}_{I0})\right){\mathit{n}}_i\otimes {\mathit{n}}_i-\frac{I_3{\lambda }_i^{-2}}{D_i}tr({\mathit{M}}_i{\mathit{A}}_{I0})\mathbf{1}\right\}\right\}\}\hfill \end{array}$$

$$\begin{array}{c}+2J^{-1}\sum _{I,J=1(I<J)}^{n_f}\sum _{k=1}^{n_n}\{\frac{\partial {Ŵ}_f(\mathit{X},I_{I04}^{(n_k)},I_{\mathrm{IJ}04}^{(n_k)})}{\partial I_{\mathrm{IJ}04}^{(n_k)}}\hfill \\ \{\sum _{i=1}^3{\lambda }_i^{n_k+2}\{\left(\frac{n_k+2}{2}tr({\mathit{M}}_i[{\mathit{a}}_{I0}\otimes {\mathit{a}}_{J0}])-\frac{I_3{\lambda }_i^{-2}}{D_i}tr({\mathit{C}}^{-1}[{\mathit{a}}_{I0}\otimes {\mathit{a}}_{J0}])\right){\mathit{n}}_i\otimes {\mathit{n}}_i\hfill \\ -\frac{I_3{\lambda }_i^{-2}}{D_i}tr({\mathit{M}}_i[{\mathit{a}}_{I0}\otimes {\mathit{a}}_{J0}])\mathbf{1}\}\}\}.\hfill \end{array}$$

(91)

The expressions of ${\mathit{M}}_i$ and ${\mathit{M}}_i$ -related entities in equation (91) can be referred to Appendix 3.

4.3. Example: the contact model of CMFs network as a fibre reinforced material

Based on the theoretical analysis and data fitting reported by Huang et al. [20,21], a neo-Hookean model of matrix is taken as the example to demonstrate the present method. In this model, the free energy functions, $Ŵ(\mathit{X},{\lambda }_1,{\lambda }_2,{\lambda }_3)$ in equation (45) (or ${Ŵ}_1(\mathit{X},{\lambda }_1,{\lambda }_2,{\lambda }_3)$ in equation (46)) and ${Ŵ}_f(\mathit{X},I_{I04}^{(n_k)},I_{\mathrm{IJ}04}^{(n_k)})$ in equation (48), are specified by the setting that

$$n_m=1,n_n=1,$$

and

$$m_1=2n_1=2.$$

And let us consider the two-families of fibres, i.e.,

$$n_f=2,$$

for specifically modelling the contact model of CMFs network.

4.3.1. Cauchy stress

Based on the above setting of the material model, the total Cauchy stress in equation (91) may be further specified as

$$\mathit{\sigma }=J^{-1}\sum _{i=1}^3{\zeta }_i{\mathit{n}}_i\otimes {\mathit{n}}_i+J^{-1}\xi \mathbf{1},$$

(92)

where the scalar functions ${\zeta }_i$ and ξ are

$$\begin{array}{c}{\zeta }_i={\beta }_i+2\sum _{I=1}^2\left\{{\overline{c}}_I^{(1)}{\lambda }_i^4\left(I_{I04}^{(2)}-1\right)\left(2tr({\mathit{M}}_i{\mathit{A}}_{I0})-\frac{I_3{\lambda }_i^{-2}}{D_i}tr({\mathit{C}}^{-1}{\mathit{A}}_{I0})\right)\right\}\hfill \\ +2{\tilde{c}}_{12}^{(1)}{\lambda }_i^4\left(I_{1204}^{(2)}-\cos ({\theta }_1-{\theta }_2)\right)\left(2tr({\mathit{M}}_i[{\mathit{a}}_{10}\otimes {\mathit{a}}_{20}])-\frac{I_3{\lambda }_i^{-2}}{D_i}tr({\mathit{C}}^{-1}[{\mathit{a}}_{10}\otimes {\mathit{a}}_{20}])\right),\hfill \end{array}$$

(93)

and

$$\begin{array}{c}\xi =-2\sum _{i=1}^3\sum _{I=1}^2\left\{{\overline{c}}_I^{(1)}\frac{I_3{\lambda }_i^2}{D_i}\left(I_{I04}^{(2)}-1\right)tr({\mathit{M}}_i{\mathit{A}}_{I0})\right\}\hfill \\ -2\sum _{i=1}^3\left\{{\tilde{c}}_{12}^{(1)}\frac{I_3{\lambda }_i^2}{D_i}\left(I_{1204}^{(2)}-\cos ({\theta }_1-{\theta }_2)\right)tr\left({\mathit{M}}_i[{\mathit{a}}_{10}\otimes {\mathit{a}}_{20}]\right)\right\},\hfill \end{array}$$

in which,

$${\beta }_i=c_1\left({\lambda }_i^2-J^{-2n}\right),$$

(94)

and from equations (58), (50) and (175) (for ${\mathit{A}}_{I0}$ and ${\mathit{a}}_{10}\otimes {\mathit{a}}_{20}$ respectively) we have

$$\begin{array}{c}{I}_{I04}^{(2)}=\frac{{\gamma }^2{\lambda }_1^2}{D_1}{(\cos {\theta }_I+{\lambda }_1\sin {\theta }_I)}^2+\frac{{\gamma }^2{\lambda }_2^2}{D_2}{(\cos {\theta }_I-{\lambda }_2\sin {\theta }_I)}^2,\hfill \\ I_{1204}^{(2)}=\frac{{\gamma }^2{\lambda }_1^2}{D_1}(\cos {\theta }_1+{\lambda }_1\sin {\theta }_1)(\cos {\theta }_2+{\lambda }_1\sin {\theta }_2)+\frac{{\gamma }^2{\lambda }_2^2}{D_2}(\cos {\theta }_1-{\lambda }_2\sin {\theta }_1)(\cos {\theta }_2-{\lambda }_2\sin {\theta }_2).\hfill \end{array}$$

Furthermore, by using the relations (175) and (177) (for ${\mathit{A}}_{I0}$ where $J=I$ and ${\mathit{a}}_{10}\otimes {\mathit{a}}_{20}$ where $I=1,J=2$ , respectively), ${\zeta }_i$ in equation (93) may be expressed explicitly as the function of ${\theta }_I$ ( $I=1,2$ ):

$$\begin{array}{c}{\zeta }_i={\beta }_i+2\sum _{I=1}^2\left\{{\overline{c}}_I^{(1)}\left(I_{I04}^{(2)}-1\right)\left({\alpha }_{i1}{\sin }^2{\theta }_I+{\alpha }_{i2}{\cos }^2{\theta }_I+{\alpha }_{i3}\sin 2{\theta }_I-1\right)\right\}\hfill \\ +2\{{\tilde{c}}_{12}^{(1)}\left(I_{1204}^{(2)}-\cos ({\theta }_1-{\theta }_2)\right)\times \hfill \\ \left({\alpha }_{i1}\sin {\theta }_1\sin {\theta }_2+{\alpha }_{i2}\cos {\theta }_1\cos {\theta }_2+{\alpha }_{i3}\sin ({\theta }_1+{\theta }_2)-\cos ({\theta }_1-{\theta }_2)\right)\},\hfill \end{array}$$

(95)

where the coefficients ${\alpha }_{\mathrm{ij}}$ ( $i,j=1,2,3$ ) are

$$\left[\begin{array}{ccc}\hfill {\alpha }_{11}\hfill & \hfill {\alpha }_{12}\hfill & \hfill {\alpha }_{13}\hfill \\ \hfill {\alpha }_{21}\hfill & \hfill {\alpha }_{22}\hfill & \hfill {\alpha }_{23}\hfill \\ \hfill {\alpha }_{31}\hfill & \hfill {\alpha }_{32}\hfill & \hfill {\alpha }_{33}\hfill \end{array}\right]=\left[\begin{array}{ccc}\hfill \frac{2{\gamma }^2{\lambda }_1^4}{D_1}\hfill & \hfill \frac{{\gamma }^2{\lambda }_1^2}{D_1}\hfill & \hfill \frac{\gamma {\lambda }_1^2}{D_1}(2\gamma {\lambda }_1+1)\hfill \\ \hfill \frac{2{\gamma }^2{\lambda }_2^4}{D_2}\hfill & \hfill \frac{{\gamma }^2{\lambda }_2^2}{D_2}\hfill & \hfill \frac{\gamma {\lambda }_2^2}{D_2}(-2\gamma {\lambda }_2+1)\hfill \\ \hfill 0\hfill & \hfill -\frac{I_3{\lambda }_3^2}{D_3}{\gamma }^2\hfill & \hfill \frac{I_3{\lambda }_3^2}{D_3}\hfill \end{array}\right].$$

4.3.2. von Mises stress as the function of $\gamma$

As aforementioned, according to growth models in Boudaoud [13], Dumais et al. [14], von Mises stress may be considered as the indication of growth criterion and measure. By using the formulation of principal stresses, von Mises stress, denoted as ${\sigma }_v$ , is specified by using equation (92) as

$${\sigma }_v=J^{-1}\sqrt{\frac{1}{2}[{({\zeta }_1-{\zeta }_2)}^2+{({\zeta }_2-{\zeta }_3)}^2+{({\zeta }_3-{\zeta }_1)}^2]}.$$

(96)

From equations (94) and (95), it is clear that, for a pair of $i,j$ $(\mathrm{\forall }i,j=1,2,3)$ , we have

$${\zeta }_i-{\zeta }_j=c_1\left({\lambda }_i^2-{\lambda }_j^2\right)$$

$$\begin{array}{c}+2\sum _{I=1}^2\left\{{\overline{c}}_I^{(1)}\left(I_{I04}^{(2)}-1\right)\left(({\alpha }_{i1}-{\alpha }_{j1}){\sin }^2{\theta }_I+({\alpha }_{i2}-{\alpha }_{j2}){\cos }^2{\theta }_I+({\alpha }_{i3}-{\alpha }_{j3})\sin 2{\theta }_I-1\right)\right\}\hfill \\ +2{\tilde{c}}_{12}^{(1)}\left(I_{1204}^{(2)}-\cos ({\theta }_1-{\theta }_2)\right)(({\alpha }_{i1}-{\alpha }_{j1})\sin {\theta }_1\sin {\theta }_2+({\alpha }_{i2}-{\alpha }_{j2})\cos {\theta }_1\cos {\theta }_2\hfill \\ +({\alpha }_{i3}-{\alpha }_{j3})\sin ({\theta }_1+{\theta }_2)-\cos ({\theta }_1-{\theta }_2)).\hfill \end{array}$$

Let us introduce

$${\theta }_1=\overline{\theta }-\mathrm{\delta \theta },{\theta }_2=\overline{\theta }+\mathrm{\delta \theta },$$

where $\overline{\theta }$ is the averaged fibre angle of the two families of fibres while $\mathrm{\delta \theta }$ the deviation due to fibre waviness (see Figure 1(b)).Then ${\zeta }_i-{\zeta }_j$ may be re-expressed as

$${\zeta }_i-{\zeta }_j=c_1{g}_{\mathrm{ij}}^{(0)}+{\overline{c}}_1^{(1)}{g}_{\mathrm{ij}}^{(1)}+{\overline{c}}_2^{(1)}{g}_{\mathrm{ij}}^{(2)}+{\tilde{c}}_{12}^{(1)}{g}_{\mathrm{ij}}^{(12)},$$

(97)

where

$$\begin{array}{c}{g}_{\mathrm{ij}}^{(0)}={\lambda }_i^2-{\lambda }_j^2,\hfill \\ g_{\mathrm{ij}}^{(1)}=2\left(I_{104}^{(2)}-1\right)\times \hfill \\ \left(({\alpha }_{i1}-{\alpha }_{j1}){\sin }^2(\overline{\theta }-\mathrm{\delta \theta })+({\alpha }_{i2}-{\alpha }_{j2}){\cos }^2(\overline{\theta }-\mathrm{\delta \theta })+({\alpha }_{i3}-{\alpha }_{j3})\sin 2(\overline{\theta }-\mathrm{\delta \theta })-1\right),\hfill \\ g_{\mathrm{ij}}^{(2)}=2\left(I_{204}^{(2)}-1\right)\times \hfill \\ \left(({\alpha }_{i1}-{\alpha }_{j1}){\sin }^2(\overline{\theta }+\mathrm{\delta \theta })+({\alpha }_{i2}-{\alpha }_{j2}){\cos }^2(\overline{\theta }+\mathrm{\delta \theta })+({\alpha }_{i3}-{\alpha }_{j3})\sin 2(\overline{\theta }+\mathrm{\delta \theta })-1\right),\hfill \\ g_{\mathrm{ij}}^{(12)}=2\left(I_{1204}^{(2)}-\cos (2\mathrm{\delta \theta })\right)\{({\alpha }_{i1}-{\alpha }_{j1})\sin (\overline{\theta }-\mathrm{\delta \theta })\sin (\overline{\theta }+\mathrm{\delta \theta })\hfill \\ +({\alpha }_{i2}-{\alpha }_{j2})\cos (\overline{\theta }-\mathrm{\delta \theta })\cos (\overline{\theta }+\mathrm{\delta \theta })+({\alpha }_{i3}-{\alpha }_{j3})\sin (2\overline{\theta })-\cos (2\mathrm{\delta \theta })\},\hfill \end{array}$$

in which

$$\begin{array}{c}{I}_{104}^{(2)}=\frac{{\gamma }^2{\lambda }_1^2}{D_1}{(\cos (\overline{\theta }-\mathrm{\delta \theta })+{\lambda }_1\sin (\overline{\theta }-\mathrm{\delta \theta }))}^2+\frac{{\gamma }^2{\lambda }_2^2}{D_2}{(\cos (\overline{\theta }-\mathrm{\delta \theta })-{\lambda }_2\sin (\overline{\theta }-\mathrm{\delta \theta }))}^2,\hfill \\ I_{204}^{(2)}=\frac{{\gamma }^2{\lambda }_1^2}{D_1}{(\cos (\overline{\theta }+\mathrm{\delta \theta })+{\lambda }_1\sin (\overline{\theta }+\mathrm{\delta \theta }))}^2+\frac{{\gamma }^2{\lambda }_2^2}{D_2}{(\cos (\overline{\theta }+\mathrm{\delta \theta })-{\lambda }_2\sin (\overline{\theta }+\mathrm{\delta \theta }))}^2,\hfill \end{array}$$

and

$$\begin{array}{c}{I}_{1204}^{(2)}=\frac{{\gamma }^2{\lambda }_1^2}{D_1}(\cos (\overline{\theta }-\mathrm{\delta \theta })+{\lambda }_1\sin (\overline{\theta }-\mathrm{\delta \theta }))(\cos (\overline{\theta }+\mathrm{\delta \theta })+{\lambda }_1\sin (\overline{\theta }+\mathrm{\delta \theta }))\hfill \\ +\frac{{\gamma }^2{\lambda }_2^2}{D_2}(\cos (\overline{\theta }-\mathrm{\delta \theta })-{\lambda }_2\sin (\overline{\theta }-\mathrm{\delta \theta }))(\cos (\overline{\theta }+\mathrm{\delta \theta })-{\lambda }_2\sin (\overline{\theta }+\mathrm{\delta \theta })).\hfill \end{array}$$

The functions $g_{\mathrm{ij}}^{(0)}$ , $g_{\mathrm{ij}}^{(1)}$ , $g_{\mathrm{ij}}^{(2)}$ , and $g_{\mathrm{ij}}^{(12)}$ in equation (97) are the function of γ, $\overline{\theta }$ and $\mathrm{\delta \theta }$ . Hence, so is von Mises stress ${\sigma }_v$ in equation (96).

4.3.3. Numerical results

It is of interest to examine the contribution of each of the functions, $g_{\mathrm{ij}}^{(0)}$ , $g_{\mathrm{ij}}^{(1)}$ , $g_{\mathrm{ij}}^{(2)}$ , and $g_{\mathrm{ij}}^{(12)}$ , to von Mises stress ${\sigma }_v$ so as to understand the characteristics of the model and to guide either the design of experiments or the fitting of the experimental data. For this purpose, the following three case studies are conducted by setting that:

$$\{\begin{array}{cc}\hfill \mathrm{Case}\mathrm{study}1:\hfill & \hfill c_1=1,{\overline{c}}_1^{(1)}={\overline{c}}_2^{(1)}={\tilde{c}}_{12}^{(1)}=0\hfill \\ \hfill \mathrm{Case}\mathrm{study}2:\hfill & \hfill {\overline{c}}_1^{(1)}={\overline{c}}_2^{(1)}=1,c_1={\tilde{c}}_{12}^{(1)}=0\hfill \\ \hfill \mathrm{Case}\mathrm{study}3:\hfill & \hfill {\tilde{c}}_{12}^{(1)}=1,c_1={\overline{c}}_1^{(1)}={\overline{c}}_2^{(1)}=0\hfill \end{array}$$

for the expression (97). The Case studies 1, 2 and 3 characterise the model behaviours of the matrix part ( $g_{\mathrm{ij}}^{(0)}$ ), uncoupled part of the fibres $(g_{\mathrm{ij}}^{(1)}+g_{\mathrm{ij}}^{(2)})$ and coupled part of the fibres $(g_{\mathrm{ij}}^{(12)})$ , respectively. The ${\sigma }_v$ obtained from the above settings are dimensionless and only serve as the indications of the characteristics of the model.

First, we examine the model behaviour of the matrix part by using Case study 1. The change of ${\sigma }_v$ in terms of γ is computed for the simply shear in a range of $\gamma \in [\tan (0^{\mathrm{o}}),\tan (80^{\mathrm{o}})]$ . The curve of ${\sigma }_v$ vs. γ is presented in Figure 4 which only shows minor nonlinear effect as a the feature of the neo-Hookean material.

Figure 4.

The typical ${\sigma }_v-\gamma$ curve of Case study 1 to represent the behaviour of the matrix material.

Second, the Case study 2 is conducted to examine the behaviours of the uncoupled part of the fibres. The curves of ${\sigma }_v$ vs. γ are presented in Figure 5 for various values of $\overline{\theta }$ and $\mathrm{\delta \theta }$ .

Figure 5.

The typical ${\sigma }_v-\gamma$ curves of Case study 2 to represent the behaviour of the uncoupled part of fibres: (a) $\overline{\theta }=0^{\mathrm{o}}$ ; (b) $\overline{\theta }=45^{\mathrm{o}}$ ; (c) $\overline{\theta }=90^{\mathrm{o}}$ .

Figure 5(a) shows the results of the specific fibre arrangements of $\overline{\theta }=0^{\mathrm{o}}$ while $\mathrm{\delta \theta }=5^{\mathrm{o}},15^{\mathrm{o}},30^{\mathrm{o}},45^{\mathrm{o}}$ , respectively. Figure 5(b) and (c) are the results of $\overline{\theta }=45^{\mathrm{o}}$ and $\overline{\theta }=90^{\mathrm{o}}$ , respectively, with the same variation of $\mathrm{\delta \theta }$ as in Figure 5(a). It is noted that, by comparison with Figure 5(b)–(c), Figure 5(a) shows both lower magnitude and lower hardening effect (except when $\mathrm{\delta \theta }=45^{\mathrm{o}}$ ). This is consistent with the experimental observation that the CMFs on the innermost layer of cell wall prefer to align in $0^{\mathrm{o}}$ with respect to the circumferential direction of the cylindrical wall [37,38]. The innermost layer of the cell wall is the location where the newly produced materials from inside the cell are deposited to maintain the volumetric growth of the cell wall. In other words, the innermost layer of cell wall is the youngest layer of the cell wall and play the key role in growth. The lower magnitude and lower hardening effect of ${\sigma }_v-\gamma$ curves mean that the layer is relatively easier to deform and grow.

Figure 5(b) shows two features of $\overline{\theta }=45^{\mathrm{o}}$ fibre arrangement. One is a clear hardening effect as γ is increasing. The other is that the variation of $\mathrm{\delta \theta }$ does not create a major effect as it is observed that all ${\sigma }_v-\gamma$ curves from $\mathrm{\delta \theta }=5^{\mathrm{o}}$ to $\mathrm{\delta \theta }=45^{\mathrm{o}}$ are closely packed together in a narrow domain in Figure 5(b) by comparison to those of $\overline{\theta }=0^{\mathrm{o}}$ and $\overline{\theta }=90^{\mathrm{o}}$ in Figure 5(a) and (c). This stability of ${\sigma }_v-\gamma$ curve may be helpful to explain how the directions of the majority of CMFs are aligned to maintain the mechanical strength and stability of the cell wall. It was reported that CMF directions gradually transfer from $0^{\mathrm{o}}$ in the innermost (and youngest) layer toward a maximum angle ( $<90^{\mathrm{o}}$ ) in the outermost (and oldest) layer [37–39]. The CMFs aligning at or close to $45^{\mathrm{o}}$ may play the major role in the cell wall mechanical property as the analysis results indicate.

The prominent feature of the results of $\overline{\theta }=90^{\mathrm{o}}$ in Figure 5(c) is that, opposite to the results of $\overline{\theta }=0^{\mathrm{o}}$ and $\overline{\theta }=45^{\mathrm{o}}$ in Figure 5(a) and (b), ${\sigma }_v$ decreases sharply for $\gamma >2.5$ when $\mathrm{\delta \theta }$ increases. This indicates that $\overline{\theta }=90^{\mathrm{o}}$ may be an instability point of ${\sigma }_v$ , which may explain the observation that CMFs usually are not exactly aligned in $90^{\mathrm{o}}$ direction, i.e., longitudinal direction of the cylindrical wall to which the wall grows [39].

Third, the Case study 3 is studied for examining the behaviours of the coupled part of the fibres. The ${\sigma }_v-\gamma$ curves of the various values of $\overline{\theta }$ and $\mathrm{\delta \theta }$ are presented in Figure 6. On one hand, Figure 6(a) shows that, for $\overline{\theta }=0^{\mathrm{o}}$ , ${\sigma }_v{|}_{\gamma >2.0}$ increases when $\mathrm{\delta \theta }$ increases. This indicates that the coupling term is stable with respect to $\mathrm{\delta \theta }$ . Hence the waviness and contacts of CMFs may effectively enhance the stability of CMFs for $\overline{\theta }=0^{\mathrm{o}}$ fibre arrangement which, as aforementioned, is on the innermost layer of cell wall where new CMFs are created/deposited. On the other hand, for $\overline{\theta }=45^{\mathrm{o}}$ and $\overline{\theta }=90^{\mathrm{o}}$ , Figure 6(b) and (c) show the behaviour opposite to that of $\overline{\theta }=0^{\mathrm{o}}$ : ${\sigma }_v{|}_{\gamma >2.5}$ decreases when $\mathrm{\delta \theta }$ increases. This instability of ${\sigma }_v$ with respect to $\mathrm{\delta \theta }$ suggests that the waviness of CMFs may loosen the cell wall when CMFs are arranged in $\overline{\theta }=45^{\mathrm{o}}$ and $\overline{\theta }=90^{\mathrm{o}}$ directions.

Figure 6.

The typical ${\sigma }_v-\gamma$ curves of Case study 3 to represent the behaviour of the coupled part of fibres: (a) $\overline{\theta }=0^{\mathrm{o}}$ ; (b) $\overline{\theta }=45^{\mathrm{o}}$ ; (c) $\overline{\theta }=90^{\mathrm{o}}$ .

For more precisely understanding the stability/instability of ${\sigma }_v$ with respect to $\mathrm{\delta \theta }$ , the ${\sigma }_v$ vs. $\mathrm{\delta \theta }$ curves of Case study 3 are shown in Figure 7. Figure 7(a)–(d) present the ${\sigma }_v$ versus $\mathrm{\delta \theta }$ curves at four specific simple shear states, $\gamma =\tan (5^{\mathrm{o}})$ , $\tan (15^{\mathrm{o}})$ , $\tan (45^{\mathrm{o}})$ , $\tan (80^{\mathrm{o}})$ , respectively. And in each figure the ${\sigma }_v-\mathrm{\delta \theta }$ curves of the three fibre arrangement, $\overline{\theta }=0^{\mathrm{o}},45^{\mathrm{o}},90^{\mathrm{o}}$ , are presented. The ${\sigma }_v-\mathrm{\delta \theta }$ curves clearly show the instability points where

$$\frac{\mathrm{d}{\sigma }_v}{\mathrm{d}(\mathrm{\delta \theta })}=0\mathrm{while}\frac{{\mathrm{d}}^2{\sigma }_v}{\mathrm{d}{(\mathrm{\delta \theta })}^2}<0.$$

Figure 7.

The typical ${\sigma }_v-\mathrm{\delta \theta }$ curves of Case study 3 to show the instability points in the behaviour of the coupled part of fibres: (a) $\gamma =\tan (5^{\mathrm{o}})$ ; (b) $\gamma =\tan (15^{\mathrm{o}})$ ; (c) $\gamma =\tan (45^{\mathrm{o}})$ ; (d) $\gamma =\tan (80^{\mathrm{o}})$ .

Figure 7(d) is consistent with what is directly observed in Figure 6 that only the $\overline{\theta }=0^{\mathrm{o}}$ fibre arrangement of CMFs is stable with respect to the changes of $\mathrm{\delta \theta }$ when $\gamma >2.5$ . On the other hand, Figure 7(a)–(c) provide more clear insight into the range of $\gamma \le 2.5$ as for this range it is difficult to visually extract information from Figure 6.

In summary, the numerical study shows that, by adding the coupling term representing the waviness and contacts of CMFs, the stiffness and stability of cell wall material may have much more variations in the variable space of $\{\overline{\theta },\mathrm{\delta \theta }\}$ . How this may correlate to cell wall property and growth needs further investigation.

5. Moiré pattern in tensile stretch of cell wall

Uniaxial or biaxial tensile stretch testing is the important method to understand material properties. In the remaining part of this study, we present the mathematical analysis of cell wall related to tensile stretch which, to the author’s best knowledge, has not been reported in the literature.

Let us consider a macroscopic tensile deformation,

$$\mathit{F}=\mathrm{diag}[{\lambda }_1,{\lambda }_2,{\lambda }_3],$$

(98)

is applied on bilayer laminae of cell wall.

Given ${\mathit{a}}_0={\mathit{E}}_1$ , two fibre directions ${\mathit{a}}_{10}^{\{1\}}$ and ${\mathit{a}}_{10}^{\{2\}}$ in ${\mathcal{B}}_0$ are specified from equation (2),

$${\mathit{a}}_{10}^{\{1\}}=\hat{\mathit{R}}(-\mathrm{\delta \theta }){\mathit{a}}_0=\left\{\begin{array}{c}\hfill \cos \mathrm{\delta \theta }\hfill \\ \hfill -\sin \mathrm{\delta \theta }\hfill \\ \hfill 0\hfill \end{array}\right\},{\mathit{a}}_{10}^{\{2\}}=\hat{\mathit{R}}(\mathrm{\Lambda }-\mathrm{\delta \theta }){\mathit{a}}_0=\left\{\begin{array}{c}\hfill \cos (\mathrm{\Lambda }-\mathrm{\delta \theta })\hfill \\ \hfill \sin (\mathrm{\Lambda }-\mathrm{\delta \theta })\hfill \\ \hfill 0\hfill \end{array}\right\}.$$

(99)

Then counterpart in ${\mathcal{B}}_t$ , ${\mathit{a}}_1^{\{1\}}$ and ${\mathit{a}}_1^{\{2\}}$ , are

$${\mathit{a}}_1^{\{1\}}=\mathit{F}{\mathit{a}}_{10}^{\{1\}}=\left\{\begin{array}{c}\hfill {\lambda }_1\cos \mathrm{\delta \theta }\hfill \\ \hfill -{\lambda }_2\sin \mathrm{\delta \theta }\hfill \\ \hfill 0\hfill \end{array}\right\},{\mathit{a}}_1^{\{2\}}=\mathit{F}{\mathit{a}}_{10}^{\{2\}}=\left\{\begin{array}{c}\hfill {\lambda }_1\cos (\mathrm{\Lambda }-\mathrm{\delta \theta })\hfill \\ \hfill {\lambda }_2\sin (\mathrm{\Lambda }-\mathrm{\delta \theta })\hfill \\ \hfill 0\hfill \end{array}\right\}.$$

(100)

Accordingly, moiré vectors in ${\mathcal{B}}_t$ are computed from equation (10) as

$${\mathit{v}}_1=\frac{-{\lambda }_2}{2\sin (\mathrm{\delta \theta })\sqrt{{\lambda }_1{\cos }^2\mathrm{\delta \theta }+{\lambda }_2{\sin }^2\mathrm{\delta \theta }}}\left\{\begin{array}{c}\hfill 0\hfill \\ \hfill 1\hfill \\ \hfill 0\hfill \end{array}\right\},$$

(101)

$${\mathit{v}}_2=\frac{1}{2\sin (\mathrm{\delta \theta })\sqrt{{\lambda }_1{\cos }^2(\mathrm{\Lambda }-\mathrm{\delta \theta })+{\lambda }_2{\sin }^2(\mathrm{\Lambda }-\mathrm{\delta \theta })}}\left\{\begin{array}{c}\hfill {\lambda }_1\sin \mathrm{\Lambda }\hfill \\ \hfill -{\lambda }_2\cos \mathrm{\Lambda }\hfill \\ \hfill 0\hfill \end{array}\right\}.$$

(102)

Consequently, the normalised moiré vectors ${\overline{\mathit{v}}}_K$ are

$${\overline{\mathit{v}}}_1=\left\{\begin{array}{c}\hfill 0\hfill \\ \hfill 1\hfill \\ \hfill 0\hfill \end{array}\right\},{\overline{\mathit{v}}}_2=\frac{1}{\sqrt{{\lambda }_1{\sin }^2\mathrm{\Lambda }+{\lambda }_2{\cos }^2\mathrm{\Lambda }}}\left\{\begin{array}{c}\hfill {\lambda }_1\sin \mathrm{\Lambda }\hfill \\ \hfill -{\lambda }_2\cos \mathrm{\Lambda }\hfill \\ \hfill 0\hfill \end{array}\right\},$$

(103)

which are demonstrated in Figure 3. The corresponding moiré angles ω in ${\mathcal{B}}_t$ can be calculated as

$$\omega =\arccos ({\overline{\mathit{v}}}_1\cdot {\overline{\mathit{v}}}_2)=-\frac{{\lambda }_2\cos \mathrm{\Lambda }}{\sqrt{{\lambda }_1{\sin }^2\mathrm{\Lambda }+{\lambda }_2{\cos }^2\mathrm{\Lambda }}}.$$

(104)

It should be emphasised that the above results are based on the assumption that both layers are subjected to the same deformation F . The situation would be much more complicated if the deformation is different from one layer to the other (see e.g. Escudero et al. [17]).

It is noted that ${\overline{\mathit{v}}}_1$ and ${\overline{\mathit{v}}}_2$ are exactly the principal directions in which the CMFs connected by the contacts in the layer 1 and layer 2, respectively, tend to separate. In the next two sections, we discuss the separation of CMFs from a shared contact and the potential growth mode resulting from the separation.

6. CMFs’ microscopic motion under macroscopic tensile stretch

6.1. CMFs’ microscopic motion problem

Zhang et al. [4] reported the experimental observations of CMFs’ motions under elastic tensile stretch which cannot be explained by the conventional passive angular reorientation model of CMFs, as stated by those authors. One of the observations was the lateral separation of two adjacent CMFs which were aligned in the transverse direction of stretch (see Figure 2(c) of Zhang et al. [4]). The observation may be modelled as a problem shown in Figure 8. The two CMFs contacted to each other before being stretched macroscopically (see Figure 8(a)). When the stretch is applied, the two CMFs are separated from each other as indicated by the evidence of the detaching of CMFs from the initial contact point (see Figure 8(b)). This may not be explained by the conventional fibre-reinforced model using the perfectly bonding assumption.

Figure 8.

Microscopic separation and relative motion of CMFs under macroscopic tensile stretch: (a) reference configuration of two CMFs, (b) solid lines represent the deformed configuration of CMFs while the dash lines are the reference configuration of CMF $\overline{P_2{P}_1}$ .

For addressing this problem, we propose a micromechanical model of the CMFs’ motion based on the preceding discussion of the contact model of cell wall. As shown in Figure 8(a), the centroids of the two CMFs are assumed to be depicted by the arcs of two circles of radius $R_0$ with the central angle in radian, $\check{\theta }(\in \{-{\check{\theta }}_0,{\check{\theta }}_0\})$ . The chords of both arcs are aligned in the ${\mathit{e}}_2$ direction. Hence, in the microscopic local coordinate system, $\check{x}$ , based on the fixed spatial frame ${\mathit{e}}_i$ ( $\mathrm{\forall }i\in \{1,2,3\}$ ) in Figure 8(a), the reference configurations of the CMFs are

$${\mathit{r}}_0=\left\{\begin{array}{c}\hfill {\check{x}}_1^0(L)\hfill \\ \hfill {\check{x}}_2^0(L)\hfill \\ \hfill {\check{x}}_3^0(L)=0\hfill \end{array}\right\}=R_0\left\{\begin{array}{c}\hfill \pm (\cos \hat{\theta }-1)\hfill \\ \hfill \sin \hat{\theta }\hfill \\ \hfill 0\hfill \end{array}\right\},$$

(105)

where $L=R_0(\check{\theta }+{\check{\theta }}_0)(\in \{0,L_0=2R_0{\check{\theta }}_0\})$ is the arc length parameter on the reference configuration.

We restrict our attention to the situations that the cross-section of the CMF was proposed to be either rectangular or ribbon-like due to the arrangement of cellulose molecules [40–43]. And it is assumed that the principal directions of the cross section of CMFs are aligned in a proper way so that the microscopic motion and bending of CMFs only take place in the principal plane ${\mathit{e}}_1\wedge {\mathit{e}}_2$ when cell wall is subjected to a macroscopic stretch.

6.2. Macro-micro-kinematical relations

Let us consider the points, $P_I$ ( $\mathrm{\forall }I\in \{1,2,3,4,5,5^{\prime }\}$ ), of the two CMFs, where $P_I$ ( $\mathrm{\forall }I\in \{1,2,3,4\}$ ) are four end points while $P_5$ and $P_{5^{\prime }}$ the two contact points of the two CMFs, respectively (see Figure 8).

In the case that the contact points, $P_5$ and $P_{5^{\prime }}$ , does not join the two CMFs together due to wall loosening, it is assumed that only the end points $P_I$ ( $\mathrm{\forall }I\in \{1,2,3,4\}$ ) are subjected to the perfectly bonding constraint. We can examine the separation of the two CMFs by calculating the relative displacement between the two contact points, $P_5$ and $P_{5^{\prime }}$ , which may be verified by the experimental data.

Let us first consider a macroscopic tensile stretch (98). The corresponding macroscopic deformation is

$$x_i={\lambda }_i{X}_i\mathrm{\forall }i\in \{1,2,3\}.$$

(106)

The perfectly bonding constraint is imposed on the end points of CMFs, $P_I$ ( $\mathrm{\forall }I\in \{1,2,3,4\}$ ), in Figure 8 so that the microscopic motion of the points is identical to macroscopic counterpart, i.e.,

$$\check{\mathit{x}}(P_I)=\mathit{x}(P_I)\mathrm{\forall }I\in \{1,2,3,4\}.$$

(107)

Hence, the displacements of the point $P_I$ relative to the point $P_J$ ( $I\ne J$ $\mathrm{\forall }I,J\in \{1,2,3,4\}$ ) in the ${\mathit{e}}_i$ -direction ( $\mathrm{\forall }i\in \{1,2,3\}$ ) are denoted as

$${\mathrm{\Delta }}_i^{\mathrm{IJ}}={\check{x}}_i(P_I)-{\check{x}}_i(P_J)-(X_i(P_I)-X_i(P_J))=({\lambda }_i-1)(X_i(P_I)-X_i(P_J)).$$

(108)

For the specified problem under consideration, we only need to consider the two relative displacements in the principal directions, ${\mathit{e}}_1$ and ${\mathit{e}}_2$ , respectively,

$${\mathrm{\Delta }}_1={\mathrm{\Delta }}_1^{13}=({\lambda }_1-1)(X_1(P_1)-X_1(P_3),$$

(109)

$${\mathrm{\Delta }}_2={\mathrm{\Delta }}_2^{12}=({\lambda }_2-1)(X_2(P_1)-X_2(P_2).$$

(110)

${\mathrm{\Delta }}_1$ is the relative displacement between two CMFs measured by the displacements of the end points of the arcs. On the other hand, ${\mathrm{\Delta }}_2$ represents the relative displacement between the two end points of the same arc along its chord direction.

Furthermore, the perfectly bonding constraint is also imposed on the rotation of CMFs at the end points $P_I$ ( $\mathrm{\forall }I\in \{1,2,3,4\}$ ). Taking the end point $P_1$ as am example, let us assume that the macroscopic reference fibre direction at $P_1$ is

$${\mathit{a}}_0(P_1)=\left\{\begin{array}{c}\hfill \cos {\theta }_0\hfill \\ \hfill \sin {\theta }_0\hfill \\ \hfill 0\hfill \end{array}\right\},$$

(111)

where the macroscopic angle ${\theta }_0$ is identical to the microscopic one ${\check{\theta }}_0$ . Then the macroscopic current fibre subjected to the perfectly bonding constraint is

$$\mathit{a}(P_1)=\mathit{F}{\mathit{a}}_0(P_1)=\left\{\begin{array}{c}\hfill {\lambda }_1\cos {\theta }_0\hfill \\ \hfill {\lambda }_2\sin {\theta }_0\hfill \\ \hfill 0\hfill \end{array}\right\}.$$

(112)

Hence the current fibre direction at $P_1$ is

$$\frac{\mathit{a}(P_1)}{||\mathit{a}(P_1)||}=\left\{\begin{array}{c}\hfill \cos \theta (P_1)\hfill \\ \hfill \sin \theta (P_1)\hfill \\ \hfill 0\hfill \end{array}\right\},$$

(113)

where

$$\theta (P_1)=\arctan \frac{{\lambda }_2\sin {\theta }_0}{{\lambda }_1\cos {\theta }_0}.$$

(114)

Due to symmetries, we have

$$\theta (P_2)=\theta (P_3)=-\theta (P_1),\theta (P_4)=\theta (P_1).$$

(115)

This determines the microscopic fibre angles, $\check{\theta }$ , as

$$\check{\theta }(P_I)=\theta (P_I)\mathrm{\forall }I\in \{1,2,3,4\}.$$

(116)

Together, the microscopic motion of CMFs (displacements and rotations), at the end points $P_I$ ( $\mathrm{\forall }I\in \{1,2,3,4\}$ ) is fully determined from the macroscopic counterpart using the perfectly bonding constraint (see equations (107) and (116)). Taking the information as the boundary conditions at the end points, we may calculate the microscopic motion of the whole CMFs.

6.3. Equations and solution of the microscopic CMFs motion

Our purpose is to calculate the microscopically relative motion between the contact points, $P_5$ and $P_{5^{\prime }}$ . For this purpose, we consider such a two-step process: step 1, both CMFs, represented by the arcs $\overline{P_2{P}_1}$ and $\overline{P_4{P}_3}$ , are subjected to (1) the relative displacement, ${\mathrm{\Delta }}_2{\mathit{e}}_2$ , between the end points of the same arc, $P_1$ ( $P_3$ ) and $P_2$ ( $P_4$ ), and (2) the specified fibre angles $\hat{\theta }(P_I)$ at all the end points $P_I$ ( $\mathrm{\forall }I\in \{1,2,3,4\}$ ); step 2, the CMF $\overline{P_2{P}_1}$ is translated as a whole by a displacement, ${\mathrm{\Delta }}_1{\mathit{e}}_1$ . In step 1, the bending deformation of both CMFs need to be calculated while the step 2 is simply to impose a rigid body translation upon the CMF $\overline{P_2{P}_1}$ . We shall discuss the step 1 in detail and then add the effect of the step 2 for the total displacement.

Due to the symmetry of the two CMFs, only the arc $\overline{P_5{P}_1}$ is formulated. We introduce a moving coordinate frame ${\hat{\mathit{e}}}_i$ ( $\mathrm{\forall }i\in \{1,2,3\}$ ) with its original point attached to $P_5$ (see Figure 8(b)). The coordinates in this frame are ${\hat{x}}_i$ while the displacements ${û}_i$ ( $\mathrm{\forall }i\in \{1,2,3\}$ ). The microscopic fibre angle on ${\hat{\mathit{e}}}_1\wedge {\hat{\mathit{e}}}_2$ plane is $\hat{\theta }$ .

Then let us consider such a rod bending problem of the step 1: the end point $P_1$ is subjected to the boundary conditions,

$${û}_2(P_1)={\mathrm{\Delta }}_2∕2,\hat{\theta }(P_1)=\theta (P_1),$$

(117)

while the point $P_5$ has the symmetric conditions,

$${û}_1(P_5)={û}_2(P_5)=0,\hat{\theta }(P_5)=0.$$

(118)

Since $P_1$ does not move in the ${\mathit{e}}_1$ -direction in the step 1, we have the relation between the motion of $P_5$ in the fixed frame ${\mathit{e}}_i$ and that of $P_1$ in the moving frame ${\hat{\mathit{e}}}_i$ (see Figure 8(b)):

$${ǔ}_1(P_5)=-{û}_1(P_1).$$

(119)

Hence the problem of solving the displacement of $P_5$ in the fixed frame is transferred into the calculation of the displacement of $P_1$ in the moving frame.

To calculate the deformation of the CMF subjected to the above all-kinematics boundary conditions (117) and (118), we may re-expressed the above bending problem as follows: while keeping all other boundary conditions the same as in equations (117) and (118), the CMF is subjected to a unknown force, f, at the point $P_1$ in the ${\mathit{e}}_2$ -direction which leads to the known displacement ${\mathrm{\Delta }}_2∕2$ in the same direction. In addition, the length of the deformed arc $\overline{P_5{P}_1}$ , should be equal to its reference length $L_0∕2$ due to the inextensibility of CMF during bending.

The deformed local configuration of the CMF in the moving frame ${\hat{\mathit{e}}}_i$ is described as

$$r=\left\{\begin{array}{c}\hfill {\hat{x}}_1(l)\hfill \\ \hfill {\hat{x}}_2(l)\hfill \\ \hfill {\hat{x}}_3(l)=0\hfill \end{array}\right\},$$

(120)

where l is the arc length parameter of the deformed CMF. Then the tangent vector of the arc is $\mathit{t}=\frac{\mathrm{d}\mathit{r}}{\mathrm{d}l}$ . We introduce a local coordinate system along the arc using a moving frame, ${\tilde{\mathit{e}}}_i(l)$ ( $\mathrm{\forall }i\in \{1,2,3\}$ ), defined as

$${\tilde{\mathit{e}}}_1=\mathit{t},{\tilde{\mathit{e}}}_2=\frac{1}{||\mathrm{d}\mathit{t}∕\mathrm{d}l||}\frac{\mathrm{d}\mathit{t}}{\mathrm{d}l},{\tilde{\mathit{e}}}_3=\frac{1}{||\mathrm{d}\mathit{t}∕\mathrm{d}l||}\mathit{t}\times \frac{\mathrm{d}\mathit{t}}{\mathrm{d}l}.$$

(121)

Then we have $\frac{\mathrm{d}\mathit{t}}{\mathrm{d}l}=\mathrm{\Omega }\times \mathit{t}$ , where Ω is the vector of the rate of rotation of the the moving frame ${\tilde{\mathit{e}}}_i$ . A tensor calculation of $t\times (\mathbf{\Omega }\times t)$ gives the expression of Ω as

$$\mathbf{\Omega }=\mathit{t}\times \frac{\mathrm{d}\mathit{t}}{\mathrm{d}l}+\mathit{t}(\mathit{t}\cdot \mathbf{\Omega }).$$

(122)

Using

$$\mathit{t}=\frac{\mathrm{d}\mathit{r}}{\mathrm{d}l}=\left\{\begin{array}{c}\hfill \frac{\mathrm{d}{\hat{x}}_1}{\mathrm{d}l}\hfill \\ \hfill \frac{\mathrm{d}{\hat{x}}_2}{\mathrm{d}l}\hfill \\ \hfill 0\hfill \end{array}\right\}=\left\{\begin{array}{c}\hfill \cos \hat{\theta }\hfill \\ \hfill \sin \hat{\theta }\hfill \\ \hfill 0\hfill \end{array}\right\},$$

(123)

it is straightforward to calculate that

$$\mathit{t}\times \frac{\mathrm{d}\mathit{t}}{\mathrm{d}l}=-\frac{\mathrm{d}\hat{\theta }}{\mathrm{d}l}{\mathit{e}}_3.$$

(124)

The negative sign is due to the increasing direction of $\hat{\theta }$ is opposite to ${\mathit{e}}_3$ . And since $\mathit{t}\cdot \mathbf{\Omega }=0$ due to no torsion deformation of CMF, we have the ‘strain’ expression of CMF:

$$\mathbf{\Omega }=-\frac{\mathrm{d}\hat{\theta }}{\mathrm{d}l}{\mathit{e}}_3,$$

(125)

or

$${\mathrm{\Omega }}_1={\mathrm{\Omega }}_2=0,{\mathrm{\Omega }}_3=-\frac{\mathrm{d}\hat{\theta }}{\mathrm{d}l}.$$

(126)

The equilibrium equation of the CMF is

$$\frac{\mathrm{d}\mathit{M}}{\mathrm{d}l}=\mathit{F}\times \mathit{t},$$

(127)

where M and F are the moment of internal stress and internal resultant force, respectively. If we use the local coordinate system ${\tilde{\mathit{e}}}_i$ , the equilibrium equation reduces to

$$\frac{\mathrm{d}{M}_3}{\mathrm{d}l}=\mathit{F}\cdot {\tilde{\mathit{e}}}_2.$$

(128)

By introducing the constitutive law of CMFs,

$$M_3=E{I}_3\left({\mathrm{\Omega }}_3+\frac{1}{R_0}\right),$$

(129)

where E is Young’s modulus of CMF and $I_3$ is the principal moment of inertia in the ${\mathit{e}}_3$ -direction, equation (128) becomes

$$-E{I}_3\frac{{\mathrm{d}}^2\hat{\theta }}{\mathrm{d}{l}^2}=\mathit{F}\cdot {\tilde{\mathit{e}}}_2.$$

(130)

As aforementioned, the unknown force $f{\mathit{e}}_2$ is assigned at points $P_1$ to give a load condition $\mathit{F}=f{\mathit{e}}_2.$ Hence the governing equation (130) becomes

$$-E{I}_3\frac{{\mathrm{d}}^2\hat{\theta }}{\mathrm{d}{l}^2}=f{\mathit{e}}_2\cdot {\tilde{\mathit{e}}}_2=-f\sin \hat{\theta }.$$

(131)

The integration of the above equation gives

$$l=\sqrt{\frac{1}{2}E{I}_3}{\int }_{\hat{\theta }(P_5)}^{\hat{\theta }}\frac{\mathrm{d}\hat{\theta }}{\sqrt{c_1-f\cos \hat{\theta }}}$$

(132)

for $\hat{\theta }\in [\hat{\theta }(P_5)=0,\hat{\theta }(P_1)=\theta (P_1)]$ . Then, using $\mathrm{d}{\hat{x}}_1∕\mathrm{d}l=\sin \hat{\theta }$ and $\mathrm{d}{\hat{x}}_2∕\mathrm{d}l=\cos \hat{\theta }$ , the deformed configuration of the arc is described as

$${\hat{x}}_1=\sqrt{\frac{1}{2}E{I}_3(c_1-f\cos \hat{\theta })∕f^2}+c_2,$$

(133)

$${\hat{x}}_2=\sqrt{\frac{1}{2}E{I}_3}{\int }_{\hat{\theta }(P_5)}^{\hat{\theta }}\frac{\cos \hat{\theta }\mathrm{d}\hat{\theta }}{\sqrt{c_1-f\cos \hat{\theta }}}+c_3.$$

(134)

The boundary condition of ${û}_2$ in equations (118) and (117) at points $P_5$ and $P_1$ yields $c_3=0$ and, since ${\hat{x}}_2={\check{x}}_2$ ,

$${\hat{x}}_2(P_1)-{\check{x}}_2^0(P_1)=\sqrt{\frac{1}{2}E{I}_3}{\int }_{\hat{\theta }(P_5)}^{\hat{\theta }(P_1)}\frac{\cos \hat{\theta }\mathrm{d}\hat{\theta }}{\sqrt{c_1-f\cos \hat{\theta }}}-\frac{H_0}{2}=\frac{{\mathrm{\Delta }}_2}{2},$$

(135)

where ${\mathit{H}}_0={\check{x}}_2^0(P_1)-{\check{x}}_2^0(P_2)=2({\check{x}}_2^0(P_1)-{\check{x}}_2^0(P_5))$ is the reference length of the chord of the arc $\overline{P_2{P}_1}$ . In addition, the in-extensibility condition of the arc $\overline{P_5{P}_1}$ reads

$$l_{\overline{P_5{P}_1}}=\sqrt{\frac{1}{2}E{I}_3}{\int }_{\hat{\theta }(P_5)}^{\hat{\theta }(P_1)}\frac{\mathrm{d}\hat{\theta }}{\sqrt{c_1-f\cos \hat{\theta }}}=\frac{L_0}{2}.$$

(136)

The two equations (135) and (136) can determine $c_1$ and f. Then, using the boundary condition of ${û}_1$ in equation (118) at point $P_5$ , we have

$${û}_1(P_5)={\hat{x}}_1(P_5)-{\check{x}}_1^0(P_5)=\sqrt{\frac{1}{2}E{I}_3(c_1-f\cos \hat{\theta }(P_5))∕f^2}+c_2-{\check{x}}_1^0(P_5)=0.$$

(137)

Here we use the relation that both r and ${\mathit{r}}_0$ of CMF₁ place $P_5$ at the original points of their coordinate systems, ${\hat{\mathit{e}}}_i$ and ${\mathit{e}}_i$ , respectively (see Figure 8(b)). Substituting the values of $c_1$ and f into (137) yields the value of $c_2$ .

Consequently, the displacement of the point $P_1$ in the ${\hat{\mathit{e}}}_1$ -direction is

$${û}_1(P_1)={\hat{x}}_1(P_1)-{\check{x}}_1^0(P_1)=\sqrt{\frac{1}{2}E{I}_3(c_1-f\cos \hat{\theta }(P_1))∕f^2}+c_2-{\check{x}}_1^0(P_1),$$

(138)

which gives the motion of $P_5$ in the fixed frame ${\mathit{e}}_i$ (see equation (119)),

$${ǔ}_1(P_5)=-{û}_1(P_1).$$

Due to symmetry, the displacement of the contact point $P_{5^{\prime }}$ is ${ǔ}_1(P_{5^{\prime }})=-{ǔ}_1(P_5)$ . Hence the motion of the contact point $P_5$ of the arc $\overline{P_2{P}_1}$ relative to the contact point $P_{5^{\prime }}$ of $\overline{P_4{P}_3}$ due to deformation in step 1 is $2{ǔ}_1(P_5)$ .

By imposing the subsequent rigid body translation ${\mathrm{\Delta }}_1{\mathit{e}}_1$ on the CMF $\overline{P_2{P}_1}$ in the step 2, the total relative motion between $P_5$ and $P_{5^{\prime }}$ is

$${\check{x}}_1(P_5)-{\check{x}}_1(P_{5^{\prime }})={\mathrm{\Delta }}_1+2{ǔ}_1(P_5),$$

(139)

In the case that ${\mathrm{\Delta }}_1>0$ and ${\mathrm{\Delta }}_2<0$ , CMFs are under compression so that ${ǔ}_1(P_5)<0$ . Hence, if two CMFs are separating from each other as shown in Figure 8(b), it requires that ${\mathrm{\Delta }}_1|$.

7. Potential growth mode

The conventional kinematical growth theory has been relying on the so-called intermediate configuration, denoted $\overline{\mathcal{B}}$ here, which is the conceptually assembly of the stress-free infinitesimal pieces of material from ${\mathcal{B}}_0$ [44]. Each piece has its own growth described by a growth deformation gradient, say ${\mathit{F}}_g$ , which generally is not compatible with its neighbours’ when the pieces are assembled together. That is the subsequent elastic deformation, denoted ${\mathit{F}}_{\mathit{e}}$ , glueing the pieces in $\overline{\mathcal{B}}$ together to form a compatible current configuration ${\mathcal{B}}_t$ in the real physics space. The intermediate configuration is considered to be fictitious [45] as it is neither observable nor measurable, which leaves a gap between a growth model and the experimental observation. In this section, by using micromechanical method, we conduct the explicit analysis in the intermediate configuration of a potential growth mode to demonstrate the possible way to bridge the gap.

In the tethering model of cell wall, the growth was considered to be dominated by the matrix both physically and mathematically. The matrix is in an amorphous state. Hence, unlike plasticity theory where crystal plasticity provides micromechanics interpretation for $\overline{\mathcal{B}}$ , to the author’s best knowledge, currently the tethering model of cell wall lacked the direct micromechanics counterpart for interpretation.

In contrast to the tethering model, the contact model suggested that CMFs play the dominant role in both the mechanical and growth of the cell wall. Hence, the motion of CMFs, which is more likely to be observed than that of matrix components, may provide the experimental evidences for the micromechanics model of growth.

Based on this consideration, according to the discussion in Sections 5 and 6, we may consider mathematically a potential growth mode of the contact model of cell wall. It is assumed that, due to the loosening of one representative contact, say point $\mathcal{O}$ in Figure 9(a), the associated CMFs separate and move away from each other, which creates a specific growth mode. By using this growth mode, we demonstrate mathematically that it is possible to conduct the micromechanics analysis of growth in $\overline{\mathcal{B}}$ .

Figure 9.

Reference configuration ${\mathcal{E}}_0$ of the representative element of a cell wall lamina according to the contact model in Figure 1(b): (a) ${\mathcal{E}}_0$ is composed of 4×6 virtual blocks referred as ‘crystals’, the arrows stand for 6 CMFs travelling from bottom to top of ${\mathcal{E}}_0$ . The contact at $\mathcal{O}$ is loosened and then the two adjacent CMFs separate and move in opposite directions. (b) the figure shows the domains ${\mathcal{D}}_1$ , ${\mathcal{D}}_2$ , ${\mathcal{D}}_3$ , the ‘crystals’ ${\mathcal{E}}_c^1$ , ${\mathcal{E}}_c^2$ , Burgers circuit $\mathcal{S}(t_0)$ , and the parallel transport of the vector u along the curve $\overline{\mathcal{C}}$ .

7.1. Analysis of the growth mode in the intermediate configuration with Burgers vector

Let us consider at the reference time $t_0$ a representative element, ${\mathcal{E}}_0(\subset {\mathcal{B}}_0)$ , of a lamina of cell wall subjected to a macroscopic tensile stretch. The representative element ${\mathcal{E}}_0$ is composed of, say, 4×6 blocks referred as ‘crystals’ as shown in Figure 9. The ‘crystals’ are defined by the characteristic length of CMFs, $l_c$ , and the waviness angle $2\mathrm{\delta \theta }$ . The growth mode is created by two steps: step 1 shear of the ‘crystals’ due to the separation and motion of CMFs, and step 2 volumetric growth of the ‘crystals’.

7.1.1. CMFs’ motion and shear of ‘crystals’

As discussed in Section 6, two CMFs may detach and separate from each other due to the loosening of the contact at the point $\mathcal{O}$ between them. It is assumed that, due to the microscopic motion, the two CMFs may move to the positions shown in Figure 10 where CMFs stop moving further due to the inextensibility of CMFs during growth. Hence the intact contact at the point $\mathcal{O}$ in Figure 9 is separated into two points in Figure 10: one still denotes $\mathcal{O}$ while the other ${\mathcal{O}}^{\prime }$ . By moving other CMFs (not associated with $\mathcal{O}$ or ${\mathcal{O}}^{\prime }$ ) in sequence, which may have different CMFs network patterns as shown in Figure 10, we obtain a growth configuration of the representative element, denoted $\overline{\mathcal{E}}(t_1)(\subset \overline{\mathcal{B}}(t_1))$ for a specific time $t_1$ ( $t_1>t_0$ to differentiate from the reference time $t_0$ ) as shown in Figure 11(a).

Figure 10.

As the contact of two CMFs at $\mathcal{O}$ is loosened and separated, domains ${\mathcal{D}}_1$ and ${\mathcal{D}}_2$ in Figure 9(b) are sheared to create a empty rhombus gap, i.e., an internal incompatibility of the representative element, indicated by the vector pointing ${\mathcal{O}}^{\prime }\to \mathcal{O}$ . Three potential CMF rearrangements are demonstrated: (a) the CMFs network keeps the same regularity as in the reference configuration ${\mathcal{E}}_0$ but the two CMFs represented by the dash lines become discontinued. (b) and (c) show two cases where all CMFs are continued but the CMFs network do not keep the original regularity. Note that in all cases four solid line CMFs follow the perfectly bonding role while two dash line CMFs do not. All CMFs are assumed inextensible.

The growth deformation gradient ${\mathit{F}}_g(\mathit{X},t_1)$ may have three different forms,

$${\mathit{F}}_g^{\mathrm{-}}(t_1)=\left[\begin{array}{ccc}\hfill 1\hfill & \hfill -{\gamma }_g^{\max }\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \end{array}\right],{\mathit{F}}_g^{\mathrm{+}}(t_1)=\left[\begin{array}{ccc}\hfill 1\hfill & \hfill {\gamma }_g^{\max }\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \end{array}\right],{\mathit{F}}_g^{\mathrm{o}}(t_1)=\left[\begin{array}{ccc}\hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \end{array}\right],$$

(140)

respectively, for the different ‘crystals’ located in three different domains, ${\mathcal{D}}_1(t_0),{\mathcal{D}}_2(t_0),{\mathcal{D}}_3(t_0)(\subset {\mathcal{E}}_0)$ in Figure 9(b). Whereby ${\gamma }_g^{\max }=\tan \mathrm{\delta \theta }$ is the maximum shear of ‘crystals’ due to the inextensibility of CMFs. It is clear that all of those deformation gradients of the ‘crystals’ in equation (140) are isochoric, just similar to crystal plasticity deformation. Hence all ‘crystals’ keep the same volume but the representative element $\overline{\mathcal{E}}(t_1)$ expands by comparison to ${\mathcal{E}}_0$ . Following the assumption in the kinematical growth theory, any ‘crystal’ subjected to ${\mathit{F}}_g$ is stress free, so is the assembly of all the ‘crystals’, $\overline{\mathcal{E}}(t_1)$ .

As shown in Figure 11(a), it is possible to define the single-valued deformation, $\overline{\mathit{x}}({\mathit{F}}_g(t_1),\mathit{X}):{\mathcal{E}}_0\to \overline{\mathcal{E}}(t_1)$ , everywhere except on the segment $\mathit{X}=X_2{\mathit{E}}_2$ where $X_2\in [0,2l_c\cos \mathrm{\delta \theta }]$ . On either sides of this segment, $\overline{x}$ may be defined in, for example, two specific ‘crystals’, ${\mathcal{E}}_c^1(\subset {\mathcal{D}}_1(t_0)\subset {\mathcal{E}}_0)$ and ${\mathcal{E}}_c^2(\subset {\mathcal{D}}_2(t_0)\subset {\mathcal{E}}_0)$ (see Figure 9(b)), respectively as:

$$\overline{\mathit{x}}({\mathit{F}}_g^{\mathrm{-}}(t_1),\mathit{X})=\left\{\begin{array}{c}\hfill X_1-{\gamma }_g^{\max }{X}_2\hfill \\ \hfill X_2\hfill \\ \hfill X_3\hfill \end{array}\right\},\overline{\mathit{x}}({\mathit{F}}_g^{\mathrm{+}}(t_1),\mathit{X})=\left\{\begin{array}{c}\hfill X_1+{\gamma }_g^{\max }{X}_2\hfill \\ \hfill X_2\hfill \\ \hfill X_3\hfill \end{array}\right\},$$

(141)

Figure 11.

Burgers vectors b in: (a) shear only configuration $\overline{\mathcal{E}}(t_1)$ , and (b) configuration $\overline{\mathcal{E}}(t)$ (dash lines) with volumetric growth imposing upon $\overline{\mathcal{E}}(t_1)$ (solid lines). Noting that the inextensibility of CMFs put a constraint on the feasible deformation of ‘crystals’ where the assumption of perfectly bonding is applicable.

which give the deformed configurations of the ‘crystals’, ${\overline{\mathcal{E}}}_c^1(t_1)$ and ${\overline{\mathcal{E}}}_c^2(t_1)(\subset \overline{\mathcal{E}}(t_1))$ , as highlighted by grey colour in Figure 11(a). Defining $\overline{x}$ for all other ‘crystals’ in the similar way may totally define $\overline{\mathcal{E}}(t_1)(\subset \overline{\mathcal{B}}(t_1))$ in the real physics space as shown in Figure 11(a). Hence we may avoid the difficulty that $\overline{\mathcal{B}}$ is a totally fictitious configuration, and could conduct further analysis in $\overline{\mathcal{B}}$ .

Noting that there is a jump of $\overline{\mathit{x}}$ across the segment $\mathit{X}=X_2{\mathit{E}}_2$ ( $X_2\in [0,l_c\cos \mathrm{\delta \theta }]$ ) between ${\mathcal{E}}_c^1$ and ${\mathcal{E}}_c^2$ . This indicates the internal incompatibility of the configuration $\overline{\mathcal{E}}(t_1)(\subset \overline{\mathcal{B}}(t_1))$ .

This incompatibility is clearly depicted by the empty rhombus gap spanning between ${\mathcal{O}}^{\prime }$ and $\mathcal{O}$ in $\overline{\mathcal{E}}(t_1)$ (Figure 11(a)). To analyse this, let us draw a closed line of Burgers circuit, $\mathcal{S}(t_0)$ , cross the contact point $\mathcal{O}$ in the reference configuration ${\mathcal{E}}_0$ (Figure 9(b)). The corresponding deformed line $\mathcal{C}(t_1)$ in $\overline{\mathcal{E}}(t_1)$ is not a closed line (Figure 11(a)). Instead, it has two end points, $\mathcal{O}$ and ${\mathcal{O}}^{\prime }$ , not meeting together, hence forms a non-vanishing vector $\overrightarrow{{\mathcal{O}}^{\prime }\mathcal{O}}$ . This can be formulated by Burgers vector [29] defined for any time t,

$$\mathit{b}(t)=-{\oint }_{\mathcal{S}(t_0)}\mathrm{d}\overline{\mathit{u}}(\mathit{X},t),$$

(142)

where $\overline{\mathit{u}}(\mathit{X},t)=\overline{\mathit{x}}({\mathit{F}}_g(t),\mathit{X})-\mathit{X}$ is the displacement vector of the material points $\mathit{X}(\in \mathcal{S}(t_0))$ . From Figure 11(a), it can be calculated that the magnitude of Burgers vector at time $t_1$ , $||\mathit{b}(t_1)||$ , is equal to the length of $\overrightarrow{{\mathcal{O}}^{\prime }\mathcal{O}}$ , i.e.,

$$||\mathit{b}(t_1)||={\overline{x}}_1({\mathit{F}}_g^{\mathrm{+}}(t_1),\mathit{X}(\mathcal{O}))-{\overline{x}}_1({\mathit{F}}_g^{\mathrm{-}}(t_1),\mathit{X}({\mathcal{O}}^{\prime }))=2{\gamma }_g^{\max }{l}_c\cos \mathrm{\delta \theta }=4l_c\sin \mathrm{\delta \theta },$$

(143)

where $\mathit{X}(\mathcal{O})=\mathit{X}({\mathcal{O}}^{\prime })={\{0,l_c\cos \mathrm{\delta \theta },X_3=\mathrm{const}.\}}^{\mathrm{T}}$ .

7.1.2. Volumetric growth of ‘crystals’

Bearing in mind that ${\mathit{F}}_g^{\mathrm{-}}$ , ${\mathit{F}}_g^{\mathrm{+}}$ and ${\mathit{F}}_g^{\mathrm{o}}$ are purely isochoric deformation of ‘crystals’. Growth may have two distinctive scenarios: (1) new materials identical to those of the ‘crystals’ may fill in the empty rhombus gap while all existing ‘crystals’ unchanged, (2) ‘crystals’ may have volumetric growth imposing upon ${\mathit{F}}_g^{\mathrm{-}}$ , ${\mathit{F}}_g^{\mathrm{+}}$ and ${\mathit{F}}_g^{\mathrm{o}}$ . Either way should increase the total volume of the representative element and reduce the magnitude of Burgers vector. We only demonstrate the second scenario as discussed below.

Let us consider a volumetric growth, which is an in-plane expansion in the ${\mathit{e}}_1$ direction,

$${\mathit{F}}_g^{\mathit{v}}(t)=\mathrm{diag}[{\lambda }_g(t),1,1],$$

(144)

where $t\ge t_1$ and ${\lambda }_g(t)\ge 1$ , imposed upon the ‘crystals’ in $\overline{\mathcal{E}}(t_1)$ already sheared by equation (140). Consequently, the total growth deformation gradient ${\mathit{F}}_g(\mathit{X},t)$ of the different ‘crystals’ in three different domains ${\mathcal{D}}_1(t_0),{\mathcal{D}}_2(t_0),{\mathcal{D}}_3(t_0)(\subset {\mathcal{E}}_0)$ in Figure 9(b), respectively, becomes:

$$\begin{array}{c}{\mathit{F}}_g^{\mathit{v}}(t){\mathit{F}}_g^{\mathrm{-}}(t)=\left[\begin{array}{ccc}\hfill {\lambda }_g(t)\hfill & \hfill -{\lambda }_g(t){\gamma }_g(t)\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \end{array}\right],{\mathit{F}}_g^{\mathit{v}}(t){\mathit{F}}_g^{\mathrm{+}}(t)=\left[\begin{array}{ccc}\hfill {\lambda }_g(t)\hfill & \hfill {\lambda }_g(t){\gamma }_g(t)\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \end{array}\right],\hfill \\ {\mathit{F}}_g^{\mathit{v}}(t){\mathit{F}}_g^{\mathrm{o}}(t)=\mathrm{diag}[{\lambda }_g(t),1,1],\hfill \end{array}$$

(145)

where ${\gamma }_g(t)$ is a function of time t satisfying the condition ${\gamma }_g(t){|}_{t=t_1}={\gamma }_g^{\max }$ . This yields the other deformed configurations of the representative element denoted as $\overline{\mathcal{E}}(t)(\subset \overline{\mathcal{B}}(t))$ for $t>t_1$ (see Figure 11(b)).

Accordingly, the total deformation, $\overline{\mathit{x}}({\mathit{F}}_g(t),\mathit{X}):{\mathcal{E}}_0\to \overline{\mathcal{E}}(t)$ , of the two specific ‘crystals’ ${\mathcal{E}}_c^1$ and ${\mathcal{E}}_c^2$ becomes

$$\overline{\mathit{x}}({\mathit{F}}_g^{\mathit{v}}(t){\mathit{F}}_g^{\mathrm{-}}(t),\mathit{X})=\left\{\begin{array}{c}\hfill {\lambda }_g(t)X_1-{\lambda }_g(t){\gamma }_g(t)X_2\hfill \\ \hfill X_2\hfill \\ \hfill X_3\hfill \end{array}\right\},$$

$$\overline{\mathit{x}}({\mathit{F}}_g^{\mathit{v}}(t){\mathit{F}}_g^{\mathrm{+}}(t),\mathit{X})=\left\{\begin{array}{c}\hfill {\lambda }_g(t)X_1+{\lambda }_g(t){\gamma }_g(t)X_2\hfill \\ \hfill X_2\hfill \\ \hfill X_3\hfill \end{array}\right\},$$

(146)

which is also highlighted by grey colour in Figure 11(b) to show the increase of the volumes.

In this case, the magnitude of Burgers vector, $\mathit{b}(t)$ of $\overline{\mathcal{E}}(t)$ (see Figure 11(b)), is

$$||\mathit{b}(t)||={\overline{x}}_1({\mathit{F}}_g^{\mathit{v}}(t){\mathit{F}}_g^{\mathrm{+}}(t),\mathit{X}(\mathcal{O}))-{\overline{x}}_1({\mathit{F}}_g^{\mathit{v}}(t){\mathit{F}}_g^{\mathrm{-}}(t),\mathit{X}({\mathcal{O}}^{\prime }))=2{\lambda }_g(t){\gamma }_g(t)l_c\cos \mathrm{\delta \theta }.$$

(147)

Noting that the inextensibility of CMFs requires that the point $\mathcal{O}$ (or ${\mathcal{O}}^{\prime }$ ) cannot move further in ${\mathit{e}}_1$ (or $-{\mathit{e}}_1$ ) direction, which exactly indicates the reduce of the magnitude of Burgers vector,

$$||\mathit{b}(t)||\le ||\mathit{b}(t_1)||.$$

(148)

That is

$${\lambda }_g(t){\gamma }_g(t)\le {\gamma }_g^{\max }.$$

(149)

To achieve the above requirement, the shear ${\gamma }_g(t)$ in $\overline{\mathcal{E}}(t)$ must decrease accordingly:

$${\gamma }_g(t)\le {\gamma }_g^{\max }.$$

(150)

This may be interpreted as the relaxation of shear deformation due to volumetric growth.

7.2. Torsion of the growth mode in material space

Besides Burgers vector, the torsion in material space, denoted T , may be used to analyse the growth mode in the configuration $\overline{\mathcal{E}}(t)(\subset \overline{\mathcal{B}}(t))$ . The relevant differential geometry may be referred to Marsden and Hughes [36].

Let us consider moving a vector u at a point ${\mathit{X}}_1$ in the ‘crystal’ ${\mathcal{E}}_c^1$ along a curve, $\overline{\mathcal{C}}$ , to a point ${\mathit{X}}_2$ in the ‘crystal’ ${\mathcal{E}}_c^2$ (see Figure 9(b)) by a parallel transport in material space ${\mathcal{B}}_0$ at time t. That is, the mapped vector ${\mathit{F}}_g({\mathit{X}}_1,t)\mathit{u}({\mathit{X}}_1)$ and ${\mathit{F}}_g({\mathit{X}}_2,t)\mathit{u}({\mathit{X}}_2)$ are parallel in the Euclidean space ${ℝ}^3$ in which $\overline{\mathcal{E}}(t)$ ( and $\overline{\mathcal{B}}(t)$ ) is embedded. This requires that

$${\mathit{F}}_g({\mathit{X}}_2,t)\mathit{u}({\mathit{X}}_2)=({\mathit{F}}_g({\mathit{X}}_1,t)+\mathrm{d}{\mathit{F}}_g({\mathit{X}}_1,t))(\mathit{u}({\mathit{X}}_1)+\mathrm{d}\mathit{u}({\mathit{X}}_1))={\mathit{F}}_g({\mathit{X}}_1,t)\mathit{u}({\mathit{X}}_1),$$

(151)

which gives the linear relation if omitting the higher-order term,

$$\mathrm{d}\mathit{u}({\mathit{X}}_1)=-{\mathit{F}}_g{({\mathit{X}}_1,t)}^{-1}\mathrm{d}{\mathit{F}}_g({\mathit{X}}_1,t)\mathit{u}({\mathit{X}}_1),$$

(152)

where $\mathrm{d}{\mathit{F}}_g({\mathit{X}}_1,t)$ takes a discrete form:

$$\mathrm{d}{\mathit{F}}_g({\mathit{X}}_1,t)={\mathit{F}}_g({\mathit{X}}_2,t)-{\mathit{F}}_g({\mathit{X}}_1,t)={\mathit{F}}_g^{\mathit{v}}(t){\mathit{F}}_g^{\mathrm{+}}(t)-{\mathit{F}}_g^{\mathit{v}}(t){\mathit{F}}_g^{\mathrm{-}}(t)=\left[\begin{array}{ccc}\hfill 0\hfill & \hfill 2{\lambda }_g(t){\gamma }_g(t)\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right].$$

(153)

Hence,

$$\mathrm{d}\mathit{u}({\mathit{X}}_1)=-\left[\begin{array}{ccc}\hfill 0\hfill & \hfill 2{\gamma }_g(t)\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right]\mathit{u}({\mathit{X}}_1).$$

(154)

On the other hand, according to the definition of connection in non-metric space, we may introduce

$$\mathrm{d}{u}_i={\mathrm{\Gamma }}_{\mathrm{kj}}^i\mathrm{d}{X}_k{u}_j,$$

(155)

where ${\mathrm{\Gamma }}_{\mathrm{kj}}^i$ is the connection of parallel transport [36], and $\mathrm{d}\mathit{X}={\mathit{X}}_2-{\mathit{X}}_1$ .

By comparison between equations (154) and (155), the only non-zero term in equation (154) indicates that

$${\mathrm{\Gamma }}_{k2}^1\mathrm{d}{X}_k=-2{\gamma }_g(t),$$

(156)

while all other components of ${\mathrm{\Gamma }}_{\mathrm{kj}}^i\mathrm{d}{X}_k$ , including ${\mathrm{\Gamma }}_{k1}^1\mathrm{d}{X}_k$ , vanish. Since $\mathrm{d}{X}_k$ is arbitrary, we have

$${\mathrm{\Gamma }}_{12}^1=c^{\mathrm{\Gamma }}{\gamma }_g(t),{\mathrm{\Gamma }}_{21}^1=0,$$

(157)

where $c^{\mathrm{\Gamma }}=-2∕l_c\sin \mathrm{\delta \theta }$ is a coefficient obtained by setting $\mathrm{d}{X}_1=l_c\sin \mathrm{\delta \theta }$ , i.e., the distance between the geometric centres of ${\mathcal{E}}_c^1$ and ${\mathcal{E}}_c^2$ . Hence, one specific component of torsion T , i.e.,

$${\mathrm{T}}_{12}^1={\mathrm{\Gamma }}_{12}^1-{\mathrm{\Gamma }}_{21}^1=c^{\mathrm{\Gamma }}{\gamma }_g(t),$$

(158)

does not vanish unless ${\gamma }_g(t)=0$ . Hence torsion T may reflect and measure the non-integrability (incompatibility) of $\overline{\mathcal{E}}(t)$ ( and $\overline{\mathcal{B}}(t)$ ). It is consistent with the analysis of Burgers vector (see equation (150)) that the incompatibility is controlled by shear ${\gamma }_g(t)$ and is not explicitly controlled by volumetric growth ${\lambda }_g(t)$ .

7.3. Discussion

Hitherto, we give the detailed analysis about the scenario in which ‘crystals’ in the representative element ${\mathcal{E}}_0(\subset {\mathcal{B}}_0)$ may have volumetric growth. It is shown that in this scenario the volume of the rhombus gap between ${\mathcal{O}}^{\prime }$ and $\mathcal{O}$ may be reduced by the volumetric growth of ‘crystals’ (see Figure 11(b)). This serves the same purpose as the other scenario where ‘crystals’ are isochoric while new materials may fill in the rhombus gap directly. Therefore, in both scenarios the rhombus gap, which is a form of incompatibility inside the representative element $\overline{\mathcal{E}}(t)$ (intra-element) created by the motion of CMFs, may reduce its size while the the representative element is expanding. On the other hand, the representative element $\overline{\mathcal{E}}(t)$ itself may have incompatibility with its neighbouring elements (inter-element), which may lead to elastic deformation to create a compatible configuration ${\mathcal{B}}_t$ .

8. Conclusion

The recent development of the research on the plant cell wall has greatly enriched the understanding of (1) the structures and roles of the components of the cell wall (e.g. CMFs, hemicellulose, pectins) and (2) the interaction between those components. Especially, the knowledge about the construction of the CMFs network and its relations with other component, which is the key to understand both the mechanical behaviours and growth of the cell wall, has shown substantial improvement. This study is aiming to propose the continuum model reflecting those recent progresses and to provide a unified mathematical and numerical framework for quantitative modelling. The proposed fibre-reinforced material model may cover both the more conventional tethering model and the recently proposed contact model of the CMFs network in a unified framework. This may provide a method for the (1) quantitative modelling and (2) model/hypothesis validation and comparison if the experimental data are available.

In addition to the anisotropy of CMFs (physical fibres), the potential anisotropies arising from the moiré pattern and the contacts between CMFs are taken into account. The moiré pattern may shed new light on the modelling of the interaction between neighbouring layers of cell wall laminae. On the other hand, the anisotropy due to the contacts may enrich the model by taking into account the shape of and the material arrangement in the single contact and the spatial arrangement of a group of contacts. Together, it is demonstrated that the extended description of anisotropies may lead to the more sophisticated cell wall model to capture the material behaviour beyond the more conventional fibre-reinforced model of cell wall.

The simple shear problem is taken as a case study to demonstrate the proposed modelling framework since shear deformation has been widely considered as one of the key mechanisms for cell wall growth. The numerical results show the consistence with the experimental observations of the CMFs arrangement in the growing plant cell wall, and also indicate the potential mechanisms which the cell wall may use to regulate its stiffness and stability for achieving a sustainable growth.

Then the tensile stretch is studied for a variety of problems. First, the moiré pattern under the tensile stretch is presented, which may indicate the principal directions to which CMFs are separated. Second, the motion of the adjacent CMFs after separating from the shared contact is discussed by using the rod theory. Third, a potential growth mode is suggested based on the separation and motion of CMFs. The two mathematical tools, Burgers vector and torsion in material space, are applied to analyse the potential incompatible growth configuration. Again, the analysis of cell wall under tensile stretch provides some insights which may be out of scope of the more conventional fibre-reinforced model of cell wall.

As the details of stresses and tangent moduli are provided, this study also provides the formulations readily to be implemented in numerical methods like finite-element method for broader applications.