Debonding of cellular structures with fibre-reinforced cell walls under shear deformation

Many natural structures are cellular solids at millimetre scale and fibre-reinforced composites at micrometre scale. For these structures, mechanical properties are associated with cell strength, and phenomena such as cell separation through debonding of the middle lamella in cell walls are key in explaining some important characteristics or behaviour. To explore such phenomena, we model cellular structures with non-linear hyperelastic cell walls under large shear deformations, and incorporate cell wall material anisotropy and unilateral contact between neighbouring cells in our models. Analytically, we show that, for two cuboid walls in unilateral contact and subject to generalised shear, gaps can appear at the interface between the deforming walls. Numerically, when finite element models of periodic structures with hexagonal cells are sheared, significant cell separation is captured diagonally across the structure. Our analysis further reveals that separation is less likely between cells with high internal cell pressure (e.g. in fresh and growing fruit and vegetables) than between cells where the internal pressure is low (e.g. in cooked or ageing plants).


Introduction
Most cellular solids are anisotropic due to the structural distribution of the cells as well as the aeolotropic properties of cell wall material. For example, at millimetre scale, wood is a cellular structure, which can be modelled as a hexagonal prismatic honeycomb, while at micrometre scale the cell walls are fibre-reinforced composites, made up of fibres of crystalline cellulose embedded in an amorphous matrix of hemicellulose and lignin, with the fibre direction nearer the cell wall axis. Features such as wood density, representing the relative quantity of the cell wall in a given volume of wood tissue, and microfibril angle play important roles in the stiffness and load-bearing capacity of this complex structure, and the impact of their variations on the tree biomechanical performance is non-trivial [1][2][3][4][5]. While wood density varies significantly among wood species, the composition and strength of the cell wall are less variable, and phenomena such as cell debonding, commonly known as "cell-peeling", consisting in "the Schematic of fruit cells in mutual contact. c Schematic of debonded fruit cells pulling apart of two halves of the cell wall which debond along the central lamella" are responsible for fast and extensive crack propagation in many wood types [6].
Cells separation through debonding of the middle lamella in cell walls is also key in explaining the property or behaviour of fruit and legumes during storage or cooking, and is decisive for the quality of food products [7][8][9] (see Fig. 1). Tissue failure by cell shearing was also observed in apples when compressed tangentially [10]. Fruit tissue is a hydrostatic structure in which individual fluid-filled cells provide resistance to compressive forces, and the fluid pressure may also influence the elastic properties of the cell walls [11,12]. Physical evidence suggests that the firmness of fruit (apple, pear, tomato) decreases during pre-harvest ripening, when the cell walls of the fruit tissue become softer, and continues to decrease in post-harvest storage due to the loss of cell-to-cell contact, even though the stiffness of the cell walls increases [13]. Ripening also involves a reduction in turgor pressure, and other physiological and mechanical factors, such as changes in cell size, wall thickness, and composition, may also contribute to changes in the strength and elasticity of the cell walls. For example, during cool storage, cells from high-maturity fruit tend to lose intercellular adhesion but maintain cell wall integrity, while cells from low-maturity fruit tend to maintain relatively high cell-to-cell adhesion but the strength of the cell wall declines, so the cells are easily ruptured [14,15]. Cellular pressure also decreases after harvest, which causes cell wall relaxation and could accelerate the process of loss of adhesion. However, both the cell wall strength and the intercellular adhesion decline as fruit enters the over-ripe stage [16].
The relevant scale at which such phenomena occur, though beyond the capacity of the human eye, can be followed by mechanical analysis and mathematical models based on microstructural evidence [17][18][19]. To effectively capture cell wall debonding, mathematical models that account for the attachment between cells, which in some structures may be sufficiently weak so that the cells separate, are required for improved predictions of large distortions and failure in cellular structures. However, obtaining suitable models that are, at the same time, physically plausible, mathematically tractable, and computationally feasible raises many theoretical and numerical challenges.
In this study, we model cellular structures with non-linear hyperelastic cell walls under large shear deformations [20] and incorporate cell wall material anisotropy and unilateral contact between neighbouring cells in our structural models [21][22][23][24][25]. The theoretical and computational challenges raised by these models range from non-linear deformation of the individual elastic cell walls, to the detection of contact and openings between individual cells, to the resolution of non-linear equations for contact.
For the cell wall material, we consider one of the most common features of many cellular solids, namely transverse isotropy, whereby the material has one axis of symmetry [26][27][28]. Theoretically, we analyse the mechanical behaviour of two cuboid cell walls in unilateral, frictionless contact and subject to generalised shear deformation, and find that, if the walls are in mutual contact in the undeformed state, then gaps can form at the interface between the sheared walls. Computationally, we extend our approach to the investigation of finite element models of periodic structures with hexagonal cells and obtain that, when the structures are sheared, gaps appear between adjacent cell walls, causing extensive cell separation diagonally across the structure. To efficiently simulate these, we extend the successive deformation decomposition procedure proposed in [29], where seamless cellular structures were treated, to structures with non-penetrative intercellular contact as follows: (i) first, a continuous deformation is computed for the entire structure, as in a compact elastic solid, where only the external boundary conditions are imposed while the cells remain in mutual contact; (ii) then, for the pre-deformed structure, the microstructural properties at individual cell level, such as the unilateral contact between cell walls and the internal cell pressure, are taken into account. For computer simulations, the two-step procedure proved significantly faster and more accurate than when the external boundary conditions and contact constraints were imposed simultaneously in a single step.
Finally, we address the question as to what is the influence of the internal cell pressure on the intercellular contact where debonding is possible? To answer this question, we impose uniform normal pressure on the internal cell walls and obtain that separation is less likely between cells with high internal cell pressure than between cells where the internal pressure is low. This is in agreement with physical observations that, under applied force, tissue from high-maturity fruit (apple, pear) breaks down into small clumps of undamaged cells, while cell walls from less mature fruit, which are relatively strongly attached to each other, will rupture [15].
Although the finite shear deformation of transversely isotropic hyperelastic solids was previously analysed in the literature, a similar study of the simultaneous shear deformation of two transversely isotropic bodies under unilateral contact has not been carried out before. Our computer simulations of cellular bodies with hyperelastic cell walls in mutual non-penetrative contact under large shear deformation and the two-step strategy which we employ to solve the multi-body contact problems more efficiently are also novel. Furthermore, our theoretical and numerical results help explain the important role of internal cell pressure in some natural cellular structures where cell debonding occurs. In Sect. 2, we briefly review the formulation of the boundary value problem with unilateral contact in finite elasticity. In Sect. 3, we analyse theoretically the mechanical behaviour of two cuboids of fibre-reinforced Mooney-Rivlin material in unilateral (frictionless) contact and subject to simultaneous generalised shear deformations, and demonstrate the appearance of gaps at the potential contact zone. In Sect. 4, we present a set of computer models representing groups of hexagonal cells in mutual unilateral contact for which we explore the effects of large strain deformations assuming that the cell wall material is a fibre-reinforced Mooney-type composite.

Hyperelastic bodies in unilateral contact
We consider a system formed from two elastic bodies made from possibly different homogeneous hyperelastic materials described by the strain energy density function W i , i = 1, 2, respectively, and situated in mutual nonpenetrative contact on part of their boundary. We denote by i , i = 1, 2, the two open, bounded, connected domains, with Lipschitz continuous boundary, occupied by the two bodies, respectively (in particular, we assume that a unit normal vector n exists almost everywhere on i = ∂ i =¯ i \ i , i = 1, 2), such that 1 ∩ 2 = ∅ (see Fig. 2).
Each body is subject to a finite elastic deformation such that J = det Gradχ i > 0 on i and χ i is injective on i . The injectivity condition guarantees that interpenetration of the matter cannot occur. However, since self-contact is permitted, the above deformation does not need to be injective on¯ i . If the spatial point x = χ i (X) corresponds to the place occupied by the particle X in the deformation χ i , i = 1, 2, then, for the deformed bodies, the equilibrium state in the absence of a body load is described in terms of the Cauchy stress by the Eulerian field equation div σ (x) = 0.
(2.1) The above governing equation is completed by a constitutive law for σ , depending on material properties, and supplemented by boundary and contact conditions. Since the domains occupied by the deformed bodies are generally unknown, the above equilibrium problem is expressed as an equivalent problem in the reference configuration, with the independent variables X ∈ 1 ∪ 2 . The corresponding Lagrangian equation of non-linear elastostatics takes the form: where P is the first Piola-Kirchhoff stress tensor, such that For the two-body system, the general boundary value problem with unilateral contact conditions is to find the displacement field u(X) = x − X, for all X ∈ 1 ∪ 2 , such that the equilibrium equation (2.2) is satisfied subject to the following conditions on the relatively disjoint, open subsets of the boundary [24,25]: (2.4) • on N , the Neumann (traction) condition 5) where N is the outward unit normal vector to N , and g N d A = τ da, where τ = σ n is the surface traction measured per unit area of the deformed state; • on C , the unilateral (non-penetrative) contact conditions, in the absence of friction forces: For the contact with a rigid obstacle, η : R 3 → R is the function describing the relative distance between the body and the surface of the obstacle, N is the outward unit normal vector to this surface (oriented towards the obstacle), d ≥ 0 is the maximum relative distance, which cannot be exceeded, between potential contact points, and g ≥ 0 is the cohesion parameter. For self-contact, η = [u(X)] · N, where [u(X)] = u(X) − u(X ) denotes the jump in the displacement at the points X = X where self-contact may occur and N and N represent the outward unit normal vectors to the boundary at those points, respectively. Then (2.6) represents the unilateral contact condition setting the permitted mutual distance between potential contact points; (2.7) is the normal force condition giving the allowed normal force acting at a contact point; (2.8) is the complementarity condition that either the maximum relative distance or the maximum normal force should be attained at each point (note that, in (2.8), if the maximum relative distance is reached, then the expression in the first bracket is equal to zero, and when the maximum normal force is achieved, the second bracket is zero); and (2.9) states that, at the points of self-contact, the action and reaction principle holds. In particular, d = 0 and g = 0 correspond to the physical case of non-penetrative cohesionless contact.

Fibre-reinforced cell walls under shear deformation
The classical problem of generalised shear deformation involves finite plane deformations of a rectangular section of a material in which straight lines parallel to the X 1 -axis are displaced relative to one another in the X 2 -direction, straight lines parallel to the X 1 -axis are displaced relative to one another in the X 2 -direction, and the straight lines parallel to the X 2 -axis in the undeformed state remain straight and parallel after the deformation. When a cuboid cell wall is subject to generalised shear, the deformation takes the form: are the Cartesian coordinates for the reference (Lagrangian, material) and the deformed (Eulerian, spatial) representation, respectively, and f is a function to be determined (see Fig. 3). In particular, for the simple shear deformation, (3.1) takes the form: where k is a positive constant (see Fig. 4) [32][33][34][35]. In this case, the straight lines parallel to the X 1 -or the X 2 -axis in the undeformed state remain straight and parallel after the deformation. For a cuboid wall of incompressible Mooney-Rivlin material under deformation (3.1), Green and Adkins [36, pp. 127-129] found that the straight lines parallel to the X 1 -axis deform in the shape of a quadratic parabola. Finite element simulations of simple and generalised shear deformations of Mooney-type materials were presented in [37,38].
Here, we consider two hyperelastic bodies made from a transversely isotropic material containing one family of extensible fibres embedded in a Mooney-Rivlin material, which is described by the strain energy density function [24,32,39] W(I 1 , I 2 , 3) where I 1 , I 2 , I 3 are the principal isotropic invariants, I 4 is the anisotropic invariant, and C 1 > 0, C 2 > 0, and C 4 > 0 are constants. If the fibres are contained in the plane (X 1 , X 2 ) and oriented in the direction where ψ ∈ [0, π/2], in the reference configuration, then the stretch of the fibre λ 4 under deformation (3.1) is given by the parameter The generalised shear of anisotropic incompressible materials was analysed by Merodio et al. [40] and Destrade et al. [41].
Assuming that the bodies are initially in contact at a common interface, we wish to determine whether they will separate if they are deformed simultaneously by the generalised shear (see Fig. 5). For the deformation (3.1), the deformation gradient is equal to  1,1,1)  (0,1,1)  (1,1,1) (1,0,1) (1),1,0) (1,0,1) where f is the derivative of f with respect to X 2 . The associated left and right Cauchy-Green tensors are, respectively, and their principal invariants are For a deforming body made of a homogeneous incompressible transversely isotropic hyperelastic material described by the strain energy function W(I 1 , I 2 , I 4 ), the Cauchy stress tensor can be represented as where β 1 = 2∂W/∂ I 1 , β −1 = −2∂W/∂ I 2 , and β 4 = 2∂W/∂ I 4 are the material response coefficients, and p is a Lagrange multiplier associated with the incompressibility constraint (3.9). In particular, if the material is reinforced with fibres that are oriented in the direction (3.4), then Then, by (3.6) and (3.11), the components of the Cauchy stress (3.10) take the form: 14) where Next, by the equilibrium equation (2.1), since (3.12)-(3.16) are satisfied, it follows that Equation (3.23) shows that p is independent of x 3 , while Eq. (3.21) implies that p is a linear function of x 1 . Next, integration of (3.22) with respect to x 2 gives where a and c are undetermined constants. By (3.21) and (3.24), The constant c remains to be obtained from the contact conditions. Integrating and substituting β 1 , β −1 , and β 4 using (3.17) yields the following cubic equation in f : where b is an arbitrary constant. In particular, if the fibres align in the X 1 -direction in the reference configuration, i.e. ψ = 0, then β 4 = 0 and (3.26) reduces to a linear equation in f , as in the isotropic case.
We now focus our investigation on the special case when the fibres align in the X 2 -direction, i.e. ψ = π/2. In this case, by (3.24) and (3.17), (3.27) and the stress components (3.12)-(3.15) take the form: (3.31) Equation (3.26) can be simplified as and since C 1 > 0, C 2 > 0, and C 4 > 0, a simple inspection of this equation shows that there is only one real function f satisfying (3.32).

Possible shear deformations
Assuming that, during deformation, the points of coordinates (X 1 , 0, 0) remain fixed and those of coordinates Substituting the above in (3.32) gives b = C 4 k 3 + (C 1 + C 2 )k. (ii) If a = 0 and the fibres are much stiffer than the matrix, so that C 4 C 1 + C 2 , then (3.32) can be written as and approximated as Solving this equation in f yields Since gaps do not appear when the two bodies deform by the same deformation, we assume f 1 (X 2 ) = k X 2 as given by (3.34) for the first cube and f 2 (X 2 ) = (3k/4)X 4/3 2 as given by (3.35) for the second cube, such that f 1 (0) = f 2 (0) = 0. We wish to verify if, in this case, the contact conditions (2.6)-(2.9) are simultaneously satisfied.
The unit normal vector at any point on the curve (0, X 2 ) is • The unilateral contact condition (2.6) between the deforming cubes is (3.36) and, since for all X 2 ∈ (0, 1), (3.37) it follows that the condition (3.36) is satisfied. • The normal forces condition (2.7) for each cube is Equivalently, by (3.28)-(3.29), Since X 1 = 0 at the interface between the two bodies in the reference configuration, the condition (3.38) is equivalent to For the first cube, a = 0 and (3.39) takes the form: (3.40) where the constant c 1 remains to be determined. For the second cube, a = C 4 k 3 and (3.39) becomes (3.41) where the constant c 2 is to be determined. • At X 2 = 0, [u] · N = 0 and the action-reaction relation (2.9) implies (3.42) • By the complementarity condition (2.8), since by (3.37) the relative displacement across the interface between the two cubes is strictly greater than zero for all X 2 ∈ (0, 1), the corresponding normal contact stress must satisfy and (3.41) can be simplified as follows: Note that, by taking k → 0, the equality in (3.46) is satisfied to the first order in k, hence (3.43) holds for sufficiently small k.
Since the contact conditions (2.6)-(2.9) with d = 0 are simultaneously satisfied for both cubes, we conclude that, under the given assumptions, it is possible for gaps to appear at the common interface between the two cubes when these are sheared simultaneously.

The cohesive effect of internal cell pressure
We now turn our attention to the normal stresses (3.38) at the external faces (X 1 , Assuming that the cells are filled and there is a normal force g 0 ≤ g exerted by the cell core on the faces of the cell walls, then, at these faces, where X 1 ∈ {−1, 1}. For the first cube, the normal stresses (3.47) are equal to that at the interface between the two cubes, i.e.
For the second cube, (3.47) takes the form: (3.49) and if k → 0, then these stresses are also equal to the normal stress at the interface between the two cubes to the first order in k.
Consequently, if g 0 = g, then gaps may appear at the interface between the two cubes, else if g 0 < g, then the cubes will remain in active contact everywhere on their interface.
In particular, when g = 0, if the cells are empty, i.e. g 0 = 0, then the sheared cell walls may debond due to gaps appearing between adjacent cell walls, whereas if the cells are filled and there is internal pressure g 0 < 0 exerted on the cell walls, then these walls remain in full contact during generalised shear deformation.

Numerical examples
In this section, we first illustrate numerically the shear deformation of two hyperelastic cuboid walls in unilateral contact and then assess computationally the shearing effects in periodic cellular structures of hyperelastic material.
The numerical examples presented here were realised within the open-source Finite Elements for Biomechanics (FEBio) software environment [42], and, in particular, the unilateral contact constraints were approximated numerically using the inbuilt FEBio finite element implementation. In the computer simulations, the transverse isotropic hyperelastic material is characterised by the strain energy function [43] W(I 1 , I 2 , where C 1 = 0.1 MPa, C 2 = 0.01 MPa, C 4 = 1 MPa, and the fibre direction is parallel to the contact surface [13].

Debonding of cuboid walls
First, we consider two cuboid walls of hyperelastic material described by (4.1), occupying the domains (−1, 0) × (−1, 1) × (0, 1) and (0, 1) × (−1, 1) × (0, 1), respectively, in the undeformed state. The bodies remain in mutual frictionless non-penetrative contact while deformed simultaneously by imposing the following external boundary conditions: the lower horizontal faces are subject to a prescribed displacement in the X 1 -direction to the left and are fixed in the X 2 -direction; the upper horizontal faces are subject to a prescribed displacement in the X 1 -direction to the right and are fixed in the X 2 -direction; and all the external faces are fixed in the X 3 -direction. The undeformed and deformed bodies are presented in Fig. 6a, where gaps across their interface are captured in the deformed state. In Fig. 6b, for the two bodies, in addition to the above boundary conditions, normal pressure is also applied on the two external side surfaces, causing the bodies to remain in full (active) contact after the deformation.

Debonding of hexagonal cells
To assess the independent influence of mechanical features such as the intercellular contact and the internal cell pressure on the collective behaviour of a group of cells under large strain deformations, we model periodic, honeycomblike structures with hexagonal cells of uniform cell size, occupying a thin square domain of (dimensionless) side one in the X 1 (horizontal)-and X 2 (vertical)-directions, and 0.1 in the X 3 (out-of-plane)-direction. Throughout the structure, adjacent cells are in mutual frictionless non-penetrative contact, and each structure is deformed by imposing the following external boundary conditions: the lower horizontal faces are subject to a prescribed displacement in the X 1 -direction to the left and are fixed in the X 2 -direction; the upper horizontal faces are subject to a prescribed displacement in the X 1 -direction to the right and are fixed in the X 2 -direction; and all the external faces are fixed in the X 3 -direction. To reduce computational cost that is mostly due to the detection of contact and openings during these deformations, we extend the successive decomposition procedure (SDP) for finite deformations of continuous cellular structures introduced in [29] to elastic structures with unilateral intercellular contact, as follows: (i) first, a continuous deformation is computed for the entire structure, as in a compact elastic solid, where only the external boundary conditions are imposed while the cells remain in mutual contact; (ii) then, for the pre-deformed structure, the microstructural properties at individual cell level, such as the unilateral contact between cell walls and the internal cell pressure, are taken into account.
Formally, if x = χ (X) ∈ R 3 denotes the finite deformation from the reference state B 0 to the final configuration B, then, in step (i), we take x = χ (X) ∈ R 3 to be a continuous deformation, and in step (ii), we assume the deformation mapping x = χ (x ) = χ χ −1 (x ) from the continuously deformed state B to the final configuration B where gaps between the cell walls can occur. This procedure is illustrated schematically in Fig. 7, where F = dχ (X)/dX > 0, F = dχ (X)/dX > 0, and F = dχ (x)/dx > 0 denote the deformation gradients for the associated invertible and orientation preserving mappings, such that F = F F .
The undeformed model structures are represented in Fig. 8, where the finite element mesh is also displayed. The sheared structures without internal cell pressure are shown in Fig. 9, where extensive cell separation diagonally across the structure is captured. For the sheared structures with uniform internal cell pressure, in Fig. 10, strong intercellular cohesion without significant gaps between cells can be observed. For these numerical simulations, the two-step procedure proved significantly faster and more accurate than when the external boundary conditions and contact constraints were imposed simultaneously in a single step.   To explain the results, for the sheared structures analysed here, assuming simple shear deformation of the entire structure in step (i), this deformation is equivalent to tension in the first principal direction combined with compression in the second principal direction, where the two principal directions are situated near the two diagonals of the deforming structure, respectively [31, pp. 101-104]. Then, when the intercellular contact is taken into account in step (ii), the cells situated along the first principal direction will separate from each other, while those along the second principal direction will be mutually compressed, as depicted in Fig. 9.
However, as demonstrated theoretically, if the internal cell pressure is sufficiently large, then the resulting normal force will act as a cohesive force across the interface between adjacent sheared walls, precluding cell separation, as shown in Fig. 10.

Conclusion
Many natural structures are cellular solids at millimetre scale and fibre-reinforced composites at micrometre scale (e.g. plant stems, vegetables, fruit). For these structures, physical properties are associated with the mechanical responses of the structural elements under applied forces, and phenomena such as cell separation through debonding of the middle lamella in cell walls are key in explaining some important characteristics or behaviour.
To explore such phenomena, we model cellular structures with non-linear hyperelastic cell walls under large shear deformations and incorporate cell wall material anisotropy and unilateral contact between neighbouring cells in our models. Analytically, we show that, for two cuboid walls in unilateral contact and subject to generalised shear, gaps can appear at the interface between the deforming walls. Numerically, when finite element models of periodic structures with hexagonal cells are sheared, significant cell separation is captured diagonally across the structure. To obtain the numerical results, we employ a successive deformation decomposition technique, in which (i) first a continuous deformation is assumed throughout the structure; (ii) then, for the pre-deformed structure, the microstructural properties at individual cell level, such as the unilateral contact between cell walls and the internal cell pressure, are taken into account. The two-step procedure proves significantly faster and more accurate than a one-step approach where the external boundary conditions and contact constraints are imposed simultaneously.
Our analysis further indicates that, under large deformations, separation is less likely between cells with high internal cell pressure than between cells where the internal pressure is low. This is in agreement with the physical observations that plant tissue under high turgor pressure (e.g. in fresh and growing fruit and vegetables) failed by cell wall rupture, whereas tissue under low turgor (e.g. in cooked or ageing plants) failed by cell separation. The way cells separate or break and release their content is critical for horticultural qualities, such as fruit texture, which is of major interest to producers around the world. As markets impose increasingly stringent quality standards, there is a demand for development of models that predict changes in fruit texture so that these processes can be managed and controlled more effectively. Even though, in most cases, it is at a cellular level that the structural basis of texture is best addressed, due to the inherent complexity and diversity of cellular structures, the explicit representation of all individual cells and their contact constraints in a structure with a very large number of cells is not feasible computationally. Nonetheless, since local changes can generate changes in the overall structural properties, the microstructural model with non-penetrative intercellular contact proposed here can be incorporated in a multiple-scale approach suitable for use in large-scale finite element computations.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.