# Coupled Bend-Twist Mechanics of Biomimetic Scale Substrate

Sanjay Dharmavaram<EMAIL_ADDRESS>Hossein Ebrahimi
<EMAIL_ADDRESS>Ranajay Ghosh<EMAIL_ADDRESS>Department of
Mathematics, Bucknell University, Lewisburg, PA  Department of Mechanical and
Aerospace Engineering, University of Central Florida, Orlando, FL

###### Abstract

We develop the mechanics of combined bending and twisting deformation of a
one-dimensional filamentous structure with protruding stiff fish-scale-like
plates embedded at an angle on the surface. We develop Cosserat kinematic
formulation along with scale contact constraints. This geometrically exact
model allows us to bypass the limitations of typical finite element
computations inherent in these systems when deflections are large. The derived
structure-property relationships reveal for the first time the combined effect
of bending and twisting on a slender fish scale inspired substrate. The model
subsumes previous models on pure bending and twisting but also shows
previously unobserved phenomena that arise due to the coupled effects of these
loads. This includes a new interpretation of kinematic locking behavior,
multiple contact regimes, asymmetric sensitivities of one curvature over the
other, and sharp transitions in the nonlinear moment-curvature and torque-
twist behaviors reflecting the complex scale engagement patterns.

###### keywords:

fish scales , biomimetic , structure-property , Cosserat , architected

††journal: Journal of the Mechanics and Physics of Solids

##  1 Introduction

Fishes are synonymous with scales, Fig. 1(a), although scales are far more
versatile in nature. They cover numerous reptiles and can also be
intermittently found in mammals such as in pangolins and armadillos [1, 2].
More interestingly, there are scale-like features in the wings of butterflies,
human hair and papillae on feline tongues [3, 4] indicating the singular
importance of the scale morphology in enhancing functions. One of the major
advantages of scales architecture is that they are are generally lightweight
additions to a substrate due to low volume fraction and yet enhance stiffness,
and multifunctionality [5, 6]. Several critical properties including
protection, locomotion, camouflaging, and thermal regulation have been
attributed to scales [7, 8, 9]. Thus, they are now intensely studied as
material templates to make armors, smart skins, soft robotics and
multifunctional surfaces [10, 11, 12, 13]. Mechanically, scales give rise to
fascinating emergent behavior such as strain stiffening, evolving
directionality and anomalous frictional response [14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24]. These behaviors can potentially aid organisms in balancing
multiple complex and often contradictory functions such as locomotion with
protection, and softness with stiffness.

Such possibilities have resulted in numerous studies in the past to understand
the nature of property enhancements brought about by scales. Early research
highlighted and confirmed the outstanding behaviors of 1-dimensional beam-like
substrates that were covered uniformly with scales using a combination of
analytical and finite element (FE) models [14, 21, 25]. These works
established precise structure-property relationships in pure bending loads for
both smooth and rough sliding between scales, and for both rigid and flexible
scales [14, 15, 25, 26]. The essential characteristics of bending behavior
such as strain stiffening and locked states were found to be universally valid
even when stiff scales were not uniformly distributed (e.g. functionally
graded [18]) or loaded under non-uniform bending [17]. The locked state is
essentially a kinematic configuration beyond with scale motion would result in
interpenetration leading to rigid behavior. In reality, somewhat before
locking commences, the contact forces on the scales would be large enough to
lead to scales deformation. Thus, bending behavior would transition from
substrate deflection (soft) to scale deformation(stiff), leading to a sharp
increase in rigidity of the overall structure. Recent studies have further
confirmed locking and nonlinear strain stiffening phenomena under pure
torsion, thereby extending their validity even in twisting [21, 22]. However,
torsion structure-property calculations while underlining the universality
also highlighted the striking differences from bending. For instance, the
kinematic locking envelopes turned out to be a complex nonlinear functions of
geometry, unlike the bending case where they were linear. In addition, the
tilt angle of the scales, Fig. 1 (b-d) had a significant impact on the nature
of locking, with some angles even precluding locking behavior. When Coulomb
friction was included between the sliding scales, more differences with the
bending case were evident [15, 22]. For instance, frictional locking was
indeed common between bending and twisting, but in twisting, the frictional
locking envelopes were highly nonlinear without closed form solution unlike
the bending case. Also, the relative energy dissipated from friction during
one cycle of loading (start to lock) was found to monotonically increase with
the friction coefficient ($\mu$) for twisting case, even though for bending
case, increasing $\mu$ did not always increase dissipated work in a cycle [22,
15]. Such differences are well anticipated as the mechanical behaviors are
dependent on the interplay of structure, load and geometry with no intrinsic
guarantees of universality. Thus, it is important to investigate individual
canonical load cases carefully and rigorously before declaring either
generalities or anomalies.

|
---|---
(a) | (b)
|
(c) | (d)
Figure 1: (a) Natural fish scales under deformation mode, adapted under CC BY
2.0 [27]. (b) Fabricated biomimetic scale metamaterials under bending
deformation. (c) Fabricated biomimetic scale metamaterials under twisting
deformation. (d) Fabricated biomimetic scale metamaterials under combined
bending and twisting deformation.

Surprisingly, the mechanics of combined bending and twisting has not been
studied. There is virtually no knowledge of the behavior of slender fish scale
inspired substrates in 3-D or spatial deflections, which is significantly more
complex but of great practical utility due to the nature of real world loads
and possibility of geometrical defects in the biomimetic substrate, leading to
cross curvatures (e.g. sagging or intrinsic twist). Prior research outlined
above show that existing models may not be universally valid, scale linearly
or even hold similar functional forms across geometry-loading combinations.
Thus entirely new models are necessary to arrive at structure-property
relationships.

We address this lacuna by developing the mechanics of a frictionless
biomimetic scale-covered beam of the geometric form shown in Fig. 1 (b-d). We
use a Cosserat rod model to reveal the nature and regimes of nonlinearity, the
interplay of bending and twisting, relationships between global (substrate)
and local (scale) deformations, and locking behavior. In the process, we
obtain the structure-property relationships. We use FE simulations and prior
results in the literature to validate our model. This paper is organized as
follows: in Sec. 2, we develop the kinematics of the filament in bending and
twisting loads, and derive the moment-curvature relationships, assuming small
strains and additive strain energies. In Sec. 3, we briefly describe the
finite element model we use to validate our theory. In Sec. 4, we conclude the
paper with results and discussion.

##  2 Contact Mechanics of Scales Engagement

We develop a kinematic description of the scale-covered substrate whose
schematic is shown Fig. 2 (a). The challenge in describing kinematics of this
system is the complicated three-dimensional configurations of the scales and
their rotation with respect to the substrate. A Cosserat kinematic
description, commonly used to model the finite elasticity of rods, is ideal in
the present context. Although the structure shown in Fig. 2 has finite width,
any variation of global strains along the width (i.e., along the x-direction)
has been neglected. We therefore model the substrate as an idealized rod.
Since Cosserat kinematics incorporates a natural director basis at every point
along the rod, the orientation of the scale relative to the substrate is easy
to describe using this basis. In this way, the scales’ orientations naturally
couple to the bending strain variables of the Cosserat description. Another
advantage of the Cosserat framework is that modeling the mechanics of these
structures becomes straightforward. We add the strain energies of the
substrate and the homogenized energy for scales; see Sec. 2.4. This work will
assume the substrate’s constitutive law to be linearly elastic which assumes
small material strains in the substrates. This has been found to be sufficient
in the practical contexts of relatively thin substrates that tend to lock
before very large strains are permitted. More complicated constitutive models
can be easily incorporated [28, 29] within our framework. Other simplifying
assumptions and the general validity of the energy formulation have been
discussed elsewhere[21, 16].

###  2.1 Global Cosserat Kinematics

Consider the schematic of the system shown in Fig. 2 (a). We model the
substrate as a Cosserat rod [28, 29] whose centroidal curve in its undeformed
reference configuration is defined by $\mathbf{R}(s)$, where $s$ represents
the arc-length along the undeformed configuration. As shown in the figure, the
planar undeformed substrate points along the z-direction, i.e.,
$\mathbf{R}(s)=s\mathbf{e}_{3}$, where $\mathbf{e}_{1}$, $\mathbf{e}_{2}$, and
$\mathbf{e}_{3}$ are the standard Cartesian basis vectors along the $x$, $y$,
and $z$ directions, respectively. We identify the reference directors of the
undeformed rod ($\mathbf{D}_{i},\;i=1,2,3$) with the cartesian basis vectors,
i.e., $\mathbf{D}_{i}=\mathbf{e}_{i}$. Let $\mathbf{r}(s)$ denote the deformed
position of the centroidal curve (shown in Fig. 2 (b)), and
$\mathbf{d}_{1}(s)$, $\mathbf{d}_{2}(s)$, and $\mathbf{d}_{3}$, the
orthonormal directors moving along the deformed rod. In the Cosserat
description, the directors of the deformed configuration are related to those
of the undeformed configuration via an orthogonal (rotation) matrix
$\mathbf{Q}(s)$:

$\mathbf{d}_{i}(s)=\mathbf{Q}(s)\mathbf{e}_{i},\text{ for }i=1,2,3.$ (1)

Differentiating (1) with respect to $s$ and substituting for $\mathbf{e}_{i}$
using the same, we obtain

$\mathbf{d}_{i}^{\prime}(s)=\mathbf{K}\mathbf{d}_{i},\text{ for }i=1,2,3,$ (2)

where

$\mathbf{K}:=\mathbf{Q}^{\prime}(s)\mathbf{Q}^{T}(s)$ (3)

is a skew-symmetric matrix associated with the bending and twisting strains in
the rod. We can associate $\mathbf{K}$ (with respect to the
$\\{\mathbf{d}_{1},\mathbf{d}_{2},\mathbf{d}_{3}\\}$ basis) with an axial
vector,
$\bm{\kappa}:=\kappa_{1}\mathbf{d}_{1}+\kappa_{2}\mathbf{d}_{2}+\kappa_{3}\mathbf{d}_{3}$,
i.e.,

$\mathbf{K}=\left(\begin{array}[]{ccc}0&-\kappa_{3}&\kappa_{2}\\\
\kappa_{3}&0&-\kappa_{1}\\\ -\kappa_{2}&\kappa_{1}&0\end{array}\right),$ (4)

and write (2) as

$\mathbf{d}_{i}^{\prime}=\bm{\kappa}\times\mathbf{d}_{i}\text{ for }i=1,2,3.$
(5)

The strain variable $\kappa_{3}$ is interpreted as the twisting strain in the
rod, while $\kappa_{1}$ and $\kappa_{2}$ as the bending strains for bending
about the $x$ and $y$ axes, respectively. Note that the definition of bending
and twisting strains in the Cosserat sense described above is distinct from
the actual 3D strains in the substrate. In Cosserat kinematics, the rod is
modeled as a 1D curve, and Cosserat strains can be interpreted as the (3D)
strains averaged over the cross-section of the substrate [28]. This
approximation is justified for thin structures. In the cases of pure bending
[14] or pure twisting [21], we identify the Cosserat strains with the twist
rates and curvature of the rod, respectively. However, the matrix description
of the strains noted in (4) is appropriate in the general case considered
here. _Henceforth in this paper, we use the term ‘strains’ strictly in this
Cosserat sense_.

We assume spatially homogeneous strain in this model analogous to prior
studies [14, 25]. This makes $\mathbf{K}$ independent of $s$ and helps set up
periodicity conditions to extract structure-property relationships. For sharp
gradients in strains or functionally graded structures, either local
periodicity [18] or discrete scale by scale approach can be used [17]. The
periodicity condition allows us to study the kinematics using a representative
volume element (RVE) of a pair of scales [14, 25].

The strain variables $\kappa_{1}$, $\kappa_{2}$, and $\kappa_{3}$ are viewed
as parameters. Equation (3) can be explicitly integrated to solve for
$\mathbf{Q}(s)$ and we obtain

$\mathbf{Q}(s)=e^{s\mathbf{K}},$ (6)

where without loss of generality we have assumed that the scale at $s=0$ is
fixed and does not change its orientation as the scales deform. That is,
$\mathbf{Q}(0)=\mathbf{I}$.

|
---|---
(a) | (b)
Figure 2: (a) Schematic of flat reference configuration of the scale-covered
rod with arc length variable along the length. (b) Schematic of deformed
configuration of the same rod under coupled bend-twist load.

Besides bending and twisting, the rod would more generally experience shearing
in $x$ and $y$ directions, and stretching in the $z$-direction (we neglect
other stretching and shear modes due to the substrate’s thin cross-section).
In Cosserat kinematics, strain variables associated with the above three modes
of deformation are given by $\nu_{1},\nu_{2}$, and $\nu_{3}$, respectively.
These are related to the deformation by

$\mathbf{r}^{\prime}(s)=\sum_{i=1}^{3}\nu_{i}(s)\mathbf{d}_{i}(s).$ (7)

If we use (1) to express $\mathbf{d}_{i}$ in terms of $\mathbf{Q}(s)$ and
$\mathbf{e}_{i}$, and after integrating the resulting equation, we obtain

$\mathbf{r}(s)=\sum_{i=1}^{3}\int_{0}^{s}\nu_{i}\mathbf{Q}(\tau)\mathbf{e}_{i}\;d\tau,$
(8)

where we have used the boundary condition $\mathbf{r}(0)=\mathbf{0}$, that is,
the rod is fixed at the origin $s=0$. In this work, we assume inextensibility
and unshearability of the rod and set $\nu_{1}=\nu_{2}=0$ and $\nu_{3}=1$.
These are reasonable approximations for the current case as soft polymers tend
to be incompressible and the cross sections are not known to distort
significantly. In the FE simulations presented in Sec. 3, we approximate
inextensibility by choosing the Poisson’s ratio $\nu=0.49$. Numerical
simulations presented in Sec. 4 seem to confirm this approximation, as we
observe no appreciable length change of the rod under the combined loading.
Neither is any significant change in cross-sectional area observed. Note that
in spite of these assumptions, there is little loss of generality. If the
model presented were to be applied for design and optimization of compressible
substrates, the above assumptions could easily be relaxed. In that case, the
position of the rod must be determined using (8) with appropriate values for
parameters $\nu_{1}$, $\nu_{2}$, and $\nu_{3}$.

###  2.2 Local Kinematics of Scales

We assume that the scales are identical in shape (rectangular with exposed
length $l$ and width $2b$), equally spaced along the rod, and are much stiffer
than the underlying substrates and thus assumed to be rigid. In addition, the
assumption of periodicity of contact (i.e., the use of RVEs) is also common
[14, 21]. Non-periodic or functionally graded systems have been investigated
previously for the cases of pure bending and pure twisting and have indicated
the salient features of these systems are still preserved under periodicity
[17, 18]. In their reference configuration, the scales are oriented parallel
to each other with the orientation described in terms of the dihedral angles
$\alpha_{0}$ and $\theta_{0}$; see Fig. 2 (a). Here, $\theta_{0}$ is the angle
made by the scale with the $z$-axis, and $\alpha_{0}$ is the scale’s tilt
angle made by its base with the $x$-axis. The separation of two adjacent
scales measured along the centroidal curve of the undeformed configuration is
$d$. As the rod deforms, the scales rigidly rotate with the rod (maintaining
orientations $\theta_{0}$ and $\alpha_{0}$ with the deformed centerline) until
contact between adjacent scales occurs. When the scales are under contact,
they continue to remain rigid, but the contact induces a change in their
orientations, _viz._ , to $\theta$ and $\alpha$ due to scale rotation on the
substrate.

We now derive conditions for contact between adjacent scales. Let $s_{i}$
represent the position of the scale $i$ along the arc-length of the undeformed
rod. Let $\Omega_{i}\subset\mathbb{R}^{3}$ be the set of all points
constituting scale $i$ in the reference configuration, and
$\omega_{i}\subset\mathbb{R}^{3}$ be the set of all points on the same, but in
the deformed configuration. These are shown schematically in Figs. 2 (a) and
(b), respectively. The midpoint of the base of the scales are located at
$\mathbf{R}_{i}:=\mathbf{R}(s_{i})$ in the reference configuration and at
$\mathbf{r}_{i}:=\mathbf{r}(s_{i})$ in the deformed configuration. Note that
for prescribed strains $\kappa_{i}$ ($i=1,2,3$), $\mathbf{r}(s_{i})$ is given
by (8) where $s_{i}=i\cdot d$. It is clear from Fig. 2 that for any given
$\mathbf{X}_{i}\in\Omega_{i}$ and $\mathbf{x}_{i}\in\omega_{i}$, the vectors
$\mathbf{X}_{i}-\mathbf{R}_{i}$ and $\mathbf{x}_{i}-\mathbf{r}_{i}$ lie on
$\Omega_{i}$ and $\omega_{i}$, respectively. Since $\mathbf{N}$ is normal to
the scale $i$ (in its reference configuration), it follows that

$\mathbf{N}\cdot[\mathbf{X}_{i}-\mathbf{R}_{i}]=0.$ (9)

Since the scales rotate rigidly under deformation, we have

$\mathbf{x}_{i}-\mathbf{r}_{i}=\mathbf{Q}_{i}\left[\mathbf{X}_{i}-\mathbf{R}_{i}\right],$
(10)

where $\mathbf{Q}_{i}:=\mathbf{Q}(s_{i})=\mathbf{Q}(i\cdot d)$ is the rotation
matrix given by (6) evaluated at $s=i\cdot d$. Two adjacent scales (say, $i=0$
and $i=1$) in their deformed configurations intersect if and only if they have
at least one common point. That is, if and only if the following equation has
at least one solution:

$\mathbf{x}_{1}-\mathbf{x}_{0}=\mathbf{0}.$ (11)

Using (10), the previous equation can be equivalently written as

$\mathbf{r}_{1}+\mathbf{Q}_{1}[\mathbf{X}_{1}-\mathbf{R}_{1}]-\mathbf{X}_{0}=\mathbf{0},$
(12a) where we have used $\mathbf{Q}_{0}=\mathbf{I}$, and
$\mathbf{R}_{0}=\mathbf{r}_{0}=\mathbf{0}$ to account for the boundary
condition of scale $i=0$. Since $\mathbf{X}_{0}$ and $\mathbf{X}_{1}$ lie on
(finite) planes $\Omega_{0}$ and $\Omega_{1}$, respectively, these points must
satisfy $\mathbf{N}\cdot\mathbf{X}_{0}=0,\text{ for
}\mathbf{X}_{0}\in\Omega_{0},$ (12b)
$\mathbf{N}\cdot\Big{[}\mathbf{X}_{1}-\mathbf{R}_{1}\Big{]}=0,\text{ for
}\mathbf{X}_{1}\in\Omega_{1}.$ (12c)

Note that solutions, $\mathbf{X}_{0}$ and $\mathbf{X}_{1}$, of (12c) are the
3D coordinates of the points of intersection of the scales, pulled back to the
reference configuration. From a computational viewpoint, while it is easy to
solve this linear system of equations, verifying if the solutions lie on
_finite_ planes $\Omega_{0}$ and $\Omega_{1}$ is not as straightforward. This
is because of two factors—finiteness of $\Omega_{1}$ and $\Omega_{2}$, and
their complicated 3D orientation in $\mathbb{R}^{3}$. We deal with both of
them by employing the following change of coordinates:

$\hat{\mathbf{X}}_{i}:=\mathbf{T}[\mathbf{X}_{i}-\mathbf{R}_{i}],\;i=0,1,$
(13)

where $\hat{\mathbf{X}}_{0}$ and $\hat{\mathbf{X}}_{1}$ are rotated coordinate
variables and

$\mathbf{T}=\left(\begin{array}[]{ccc}\cos\alpha&0&\sin\alpha\\\ 0&1&0\\\
-\sin\alpha&0&\cos\alpha\end{array}\right),$ (14)

is a rotation matrix. The transformation rule (13) maps $\Omega_{i}$ to
$\hat{\Omega}_{i}:=\mathbf{T}(\Omega_{i})$, where the latter’s projection on
the $\hat{X}_{i}$-$\hat{Z}_{i}$ plane is shown in gray in Fig. (3). Thus,
$\mathbf{T}$ rotates the scales such that

$\hat{\Omega}_{0}=\hat{\Omega}_{1}=\mathcal{B}:=[-b,b]\times[0,\infty]\times[0,l\cos\theta].$
(15)

Figure 3: Schematic showing the coordinate system for $\hat{\mathbf{X}}_{i}$.
The dashed lines represent the edges of the substrate and the gray region is
the projection of the scale on the $X-Z$ plane. For each scale $i$, matrix
$\mathbf{T}$ is a coordinate transformation from $XYZ$ to
$\hat{X}_{i}\hat{Y}_{i}\hat{Z}_{i}$.

Rewriting (12c) in terms of $\hat{\mathbf{X}}_{0}$ and $\hat{\mathbf{X}}_{1}$
and using (13), we obtain the following equivalent conditions for the
intersection of the two adjacent scales:

$\mathbf{T}^{T}\hat{\mathbf{X}}_{0}-\mathbf{Q}_{1}\mathbf{T}^{T}\hat{\mathbf{X}}_{1}=\mathbf{r}_{1},$
(16a) $\mathbf{N}^{T}\mathbf{T}^{T}\hat{\mathbf{X}}_{0}=0,$ (16b)
$\mathbf{N}^{T}\mathbf{T}^{T}\hat{\mathbf{X}}_{1}=0,$ (16c)

where $\hat{\mathbf{X}}_{0},\hat{\mathbf{X}}_{1}\in\mathcal{B}$ and
superscript $T$ denotes matrix transpose. Thus, transformed variables
$\hat{\mathbf{X}}_{i}$ lie in a simple rectangular domain $\mathcal{B}$.
Transforming the variables in this manner also resolves the issue of
finiteness noted above because it is straightforward to check if solutions lie
in $\mathcal{B}$ by checking the bounds on the variables.

###  2.3 Contact Conditions

Observe that (16c) is an under-determined linear system with five equations in
six unknowns. If $\hat{\mathbf{X}}_{0}$ and $\hat{\mathbf{X}}_{1}$ remain
unrestricted (i.e., $\hat{\mathbf{X}}_{i}\in\mathbb{R}^{3}$, instead of
$\mathcal{B}$), then we have two possibilities—no solutions, when the two
planes containing the scales are parallel, or infinitely many solutions with a
line of intersection. However, since the scales constitute finite planes,
since $\hat{\mathbf{X}}_{0}$ and $\hat{\mathbf{X}}_{1}$ are restricted to
$\mathcal{B}$. When finite planes intersect, there are the following
possibilities: a) no solution, when the planes are either parallel or do not
intersect within the domain, b) infinitely many solutions when the two planes
interpenetrate each other with a line-segment intersection c) _unique or
infinitely_ many solutions (depending on the relative orientation of the
scales), but the scales _do not interpenetrate_ each other. The third case is
of interest here as they correspond to physically realizable configurations of
scales making contact. A schematic for these configurations is shown in Fig. 4
with intersections (if any) shown in red.

Figure 4: Various possibilities for two finite rectangles to intersect. The
left-most figure shows no contact between scales, the next shows
interpenetration of scales. The last three are a non-exhaustive sampling of
possible modes of contact. Intersections are shown in red.

There are various ways in which the third case of intersecting yet non-
penetrating scales may be realized. Observe that in any of these
configurations there is _at least one point of intersection_ that falls in the
following two cases: 1) it is either a corner of one of the scales (_corner
contact_) or 2) is a point of intersection of two edges of different scales
(_edge contact_). See last three figures in Fig. 4). To find the intersection
points for the first case, we prescribe corner coordinates $\hat{X}_{i}$ and
$\hat{Z}_{i}$ (for scale $i$) and for the second case, any two coordinates but
from different scales, e.g., $\hat{X}_{i}$ and $\hat{Z}_{j}$ ($i\neq j$) are
prescribed. We use the bounds for the corresponding variables appearing in
(15) to fix these values. Since $\hat{Y}_{i}$ is unbounded we do not fix this
coordinate. The two cases are not necessarily mutually exclusive, as some
configurations can be both. For instance, when corner contact happens on an
edge; see third schematic in Fig. 4. The cases described above do not
automatically exclude interpenetrating configurations. For example, the second
schematic in Fig. 4 is an edge-contact. In computations, we exclude such cases
by explicitly checking for interpenetration.

There are eight possible corner contacts depending on which one of the eight
edges (between the two scales) makes contact with the other. And there are
sixteen edge contacts depending which one of the four edges of first scale
makes contact with the four of the other. Of these twenty four possibilities
most correspond to configurations with extreme twists that lead to appreciable
material nonlinearities or even cause the scales to debond from the substrate.
Therefore, we only consider the following four cases, shown schematically in
Fig. 5:

|  |  |
---|---|---|---
(a) | (b) | (c) | (d)
Figure 5: Modes of contact for an RVE considered in this work: (a) Corner
contact (red) (b) Top-left edge contact (black) (c) Top-right edge (green) (d)
Top edge (blue). Contact point is highlighted as a red dot

  1. 1.

Corner contact: Either of the two top corners of scale ‘$0$’ makes contact
with scale ‘$1$’. For such a configuration, we set
$(\hat{X}_{0},\hat{Z}_{0})=(\pm b,l\cos\theta)$. This is shown as red in our
schematic Fig. 5 (a). In practice, $\hat{X}_{0}=-b$ is found only when
$\alpha<0$.

  2. 2.

Top-Left edge contact: This happens when the top edge of scale ‘0’ makes
contact with left edge of scale ‘1’. For this, we set
$(\hat{X}_{1},\hat{Z}_{0})=(-b,l\cos\theta)$. This is shown as black in our
schematic; see Fig. 5 (b).

  3. 3.

Top-Right edge contact: This happens when the top edge of scale ‘0’ makes
contact with left edge of scale ‘1’. For this, we set
$(\hat{X}_{1},\hat{Z}_{0})=(b,l\cos\theta)$. This is shown as green in our
schematic; see Fig. 5 (c).

  4. 4.

Top-edge contact: This happens when the top edge of scale ‘0’ makes contact
with top edge of scale ‘1’. For this, we set
$(\hat{Z}_{0},\hat{Z}_{1})=(l\cos\theta,l\cos\theta)$. This is shown as blue
in our schematic; see Fig. 5 (c).

These color codes will be used later in Sec. 4 . All the above cases can be
expressed mathematically as the conditions

$\mathbf{a}_{0}^{T}\hat{\mathbf{X}}_{0}+\mathbf{a}_{1}^{T}\hat{\mathbf{X}}_{1}=A,$
(17a)
$\mathbf{b}_{0}^{T}\hat{\mathbf{X}}_{0}+\mathbf{b}_{1}^{T}\hat{\mathbf{X}}_{1}=B,$
(17b)

where $\mathbf{a}_{0},\;\mathbf{a}_{1},\;\mathbf{b}_{0},\;\mathbf{b}_{1}$ are
3D vectors that select the appropriate coordinates of $\hat{\mathbf{X}}_{0}$
and $\hat{\mathbf{X}}_{1}$ that we wish to fix, and $A$ and $B$ are the bounds
on the coordinates. For example, for the corner contact case, i.e., Case 1,
noted above with $(\hat{X}_{0},\hat{Z}_{0})=(b,l\cos\theta)$, (17b), we have
$\mathbf{a}_{0}=(1,0,0)\;\mathbf{a}_{1}=(0,0,0)^{T}$ and
$\mathbf{b}_{0}=(0,0,1),\;\mathbf{b}_{1}=(0,0,0)^{T}$ (which pick coordinates
$\hat{X}_{0}$ and $\hat{Z}_{0}$, respectively), and $A=b$, $B=l\cos\theta$.
For Case 2, $(\hat{X}_{0},\hat{X}_{1})=(b,b)$, so
$\mathbf{a}_{0}=(1,0,0),\;\mathbf{a}_{1}=(0,0,0)^{T}$,
$\mathbf{b}_{0}=(0,0,0),\;\mathbf{b}_{1}=(1,0,0)^{T}$, $A=b$, and $B=b$.

To obtain an implicit expression for the dependence of $\mathbf{K}$ (and
therefore, $\kappa_{1}$, $\kappa_{2}$ and $\kappa_{3}$) on $\theta$, we solve
(16c) simultaneously with (17b). We do this in the following steps. First, we
solve (16a) for $\hat{\mathbf{X}}_{1}$, i.e.,
$\hat{\mathbf{X}}_{1}=\mathbf{T}\mathbf{Q}_{1}^{T}\mathbf{T}^{T}\hat{\mathbf{X}_{1}}-\mathbf{T}\mathbf{Q}_{1}^{T}\mathbf{r}_{1}$.
We then plug the previous expression for $\hat{\mathbf{X}}_{1}$ into (17a) and
(17a) to obtain equations only involving $\hat{\mathbf{X}}_{0}$. Taking the
resulting equations together with (16b), we obtain a $3\times 3$ system of
linear equations $\mathbf{A}\hat{\mathbf{X}}_{0}=\mathbf{c}$, where

$\mathbf{A}=\left(\begin{array}[]{c}(\mathbf{T}\mathbf{N})^{T}\\\
(\mathbf{a}_{0}+\mathbf{T}\mathbf{Q}_{1}\mathbf{T}^{T}\mathbf{a}_{1})^{T}\\\
(\mathbf{b}_{0}+\mathbf{T}\mathbf{Q}_{1}\mathbf{T}^{T}\mathbf{b}_{1})^{T}\end{array}\right),\;\mathbf{c}=\left(\begin{array}[]{c}0\\\
A+\mathbf{a}_{1}^{T}\mathbf{T}\mathbf{Q}_{1}^{T}\mathbf{r}_{1}\\\
B+\mathbf{b}_{1}^{T}\mathbf{T}\mathbf{Q}_{1}^{T}\mathbf{r}_{1}\\\
\end{array}\right).$ (18)

Plugging the solution $\hat{\mathbf{X}}_{0}=\mathbf{A}^{-1}\mathbf{c}$ along
with the solution
$\hat{\mathbf{X}}_{1}=\mathbf{T}\mathbf{Q}_{1}^{T}\mathbf{T}^{T}\hat{\mathbf{X}_{1}}-\mathbf{T}\mathbf{Q}_{1}^{T}\mathbf{r}_{1}$
in (16a) we obtain the following implicit dependence of $\theta$ in terms of
$\mathbf{K}$:

$f(\mathbf{K},\theta)=\mathbf{N}^{T}\mathbf{Q}_{1}^{T}\mathbf{T}^{T}\mathbf{A}^{-1}\mathbf{c}-\mathbf{N}^{T}\mathbf{Q}_{1}^{T}\mathbf{r}_{1}=0.$
(19)

The dependence of $f$ on $\mathbf{K}$ and $\theta$ can be gleaned by observing
that $\mathbf{Q}_{1}$ and $\mathbf{r}_{1}$ depend on $\mathbf{K}$ ((cf. (6)
and (8)), while $\mathbf{N}$, $A$ and $B$ depend on $\theta$. Parameters
$\alpha$ appears in $\mathbf{T}$, and $l$ and $b$ appear in $A$, and $B$.
Equation (19) is a highly nonlinear ‘bridging’ law linking the local to the
global kinematics completing the multiscale kinematics description of the
system. We solve this equation numerically to obtain the dependence of
$\theta$ on the bending strains, $\kappa_{i}$. Solutions corresponding to
interpenetrating configurations are discarded.

###  2.4 Mechanics of Biomimetic Scale Slender Substrate

In order to understand the mechanics of this structure, we take recourse to
energy balance between the global loads and local deformation. We model the
underlying substrate as an Euler elastica whose strain energy is given by:

$\mathcal{E}_{beam}=\int_{0}^{L}\Big{[}\frac{1}{2}B_{1}\kappa_{1}^{2}+\frac{1}{2}B_{2}\kappa_{2}^{2}+\frac{1}{2}D\kappa_{3}^{2}\Big{]}\;ds.$
(20)

where $B_{1},B_{2}$ and $D$ are respectively the bending and twisting
rigidities. For beams with circular cross sections,
$B_{1}=EI_{1},B_{2}=EI_{2}$ and $D=GJ$ where $E$ is the Young’s modulus,
$I_{1},I_{2}$ are respective area moments, $G$ is the shear modulus, and $J$
is the polar moment. In our case, two complications arise. The embedding of
the scales lead to an inclusion effect on the beam, resulting in additional
stiffness (composite effect). For bending deformation, we postulate a scaling
parameter $C_{f}$ to the bending rigidities, i.e.
$B_{1}=C_{f}EI_{1},B_{2}=C_{f}EI_{2}$. The values of the parameter $C_{f}$ can
be obtained using FE simulations [21]. For the twisting deformation, a non-
circular cross section also introduces warping effect. Although warping is not
substantial to affect the kinematics of the scales engagement, it does have an
effect on the energy. This is addressed using another multiplicative factor
$C_{w}$, which is readily available in handbooks [30] and often scaled with
respect to $EI_{1}$ rather than $GJ$. Thus, the twisting rigidity takes the
form $D=C_{f}C_{w}EI_{1}$.

To include the contribution from the scales, we note that a scale’s rotations
is captured by the angles $\theta$ and $\alpha$. As the substrate deforms
under applied strain, the scales rotate freely until a critical threshold of
curvatures is reached when the scales engage. Let
$\Gamma_{e}\subset\mathbb{R}^{3}$ denote the bending strains
$(\kappa_{1},\kappa_{2},\kappa_{3})$, for which the scales are engaged. This
region can be determined by solving (19) for $\kappa_{1}$, $\kappa_{2}$ and
$\kappa_{3}$ for which $\theta>\theta_{0}$ and when scales make contact. We
determined it numerically. Clearly, the region will depend on the parameters
chosen. In Fig. 10 of Sec. 4, we show a plot for a slice of $\Gamma_{e}$ along
$\kappa_{2}=0$. Since the scales are built into the substrate, under
engagement (i.e., when scales make contact), their rotations are resisted by
the substrate. We model the substrate resistance using linear torsional
springs with the elastic energy stored in the springs for each scale given by:

$\mathcal{E}_{scale}(\theta,\alpha)=\frac{1}{2d}\Big{[}K_{\theta}(\theta-\theta_{0})^{2}+K_{\alpha}(\alpha-\alpha_{0})^{2}\Big{]}H_{\Gamma_{e}}(\kappa_{1},\kappa_{2},\kappa_{3}),$
(21)

where $K_{\theta}$ and $K_{\alpha}$ are the spring constants of the torsional
springs, and $H_{\Gamma_{e}}$ is the indicator function on $\Gamma_{e}$ (i.e.,
$H_{\Gamma_{e}}(\kappa_{1},\kappa_{2},\kappa_{3})=1,\text{ if
}(\kappa_{1},\kappa_{2},\kappa_{2})\in\Gamma_{e}$, and zero, otherwise). The
total energy per RVE can be additively written as:

$\mathcal{E}(\kappa_{1},\kappa_{2},\kappa_{3})=\mathcal{E}_{beam}(\kappa_{1},\kappa_{2},\kappa_{3})+\mathcal{E}_{scale}(\theta,\alpha).$
(22)

Our FE simulations show that the change in $\alpha$ from $\alpha_{0}$ is
minimal even when the scales are engaged. This observation is in agreement
with our earlier findings for the cases of pure bending and twisting [14, 21].
We therefore drop its dependence from $\mathcal{E}$ and henceforth fix
$\alpha$ to $\alpha_{0}$. Note that in light of (19), $\theta$ is itself a
function of $(\kappa_{1},\kappa_{2},\kappa_{3})$, hence the dependence of
$\mathcal{E}$ on the same.

$K_{\theta}$ is related to the Young’s modulus of the substrate ($E$), scale
thickness ($t_{s}$), inclusion length $(L)$ and $\theta_{0}$. As we have shown
in [14, 21] the following non-dimensional scaling exists:

$\frac{K_{\theta}}{Et_{s}^{2}}=C_{B}\left(\frac{L}{t_{s}}\right)^{n}f(\theta_{0}),$
(23)

where $n$ is non-dimensionless constant that we estimate using $FE$
simulations, $C_{B}$ is a constant, and $f(\theta_{0})$ is function of angle
$\theta$. For the results presented below, these were $C_{B}=3.62$, $n=1.55$,
$f(\theta_{0})\approx 1$ [21, 22].

The moments in the three directions are computed by differentiating (22) with
respect to $\kappa_{1}$, $\kappa_{2}$, and $\kappa_{3}$:

$M_{1}=C_{f}{E}I_{1}\kappa_{1}+\frac{K_{\theta}}{d}(\theta-\theta_{0})\frac{\partial\theta}{\partial\kappa_{1}}H_{\Gamma_{e}}(\kappa_{1},\kappa_{2},\kappa_{3}),$
(24a)
$M_{2}=C_{f}{E}I_{2}\kappa_{2}+\frac{K_{\theta}}{d}(\theta-\theta_{0})\frac{\partial\theta}{\partial\kappa_{2}}H_{\Gamma_{e}}(\kappa_{1},\kappa_{2},\kappa_{3}),$
(24b)
$M_{3}=C_{f}C_{w}GI_{1}\kappa_{3}+\frac{K_{\theta}}{d}(\theta-\theta_{0})\frac{\partial\theta}{\partial\kappa_{3}}H_{\Gamma_{e}}(\kappa_{1},\kappa_{2},\kappa_{3}).$
(24c)

These relations can be computed numerically to obtain the moment-curvature
relationships. Note that the derivative of the indicator function,
$H_{\Gamma_{e}}$ is the (surface) Dirac Delta distribution which is zero
everywhere except on the boundary of $\Gamma_{e}$. This term does not appear
in (24c) because it is multiplied by $\theta-\theta_{0}$, which regularizes to
zero precisely on the boundary of $\Gamma_{e}$.

##  3 Finite Element Analysis

We build Finite Element (FE) models on ABAQUS/CAE 2017 (Dassault Systèmes) to
validate the analytical model developed above. In the simulations, we model
the substrate and the scales as 3D deformable solids. The substrate is a
rectangular prismatic member (length $L_{B}=200$ mm) onto which a row of 19
identical scales are embedded on one side. In this assembly, scales are spaced
$d=10$ mm apart, oriented with angle of $\theta_{0}=5^{\circ}$ with respect to
the substrate’s top surface, and angle of $\alpha_{0}=30^{\circ}$ with respect
to the substrate’s rectangular cross section. Since our analytical model
applies to an RVE, where strains are constant, we chose our substrate to be
sufficiently long so as to avoid edge effects. The substrate was modeled as a
linear elastic material with an elastic modulus $E_{B}=2.5$ MPa and Poisson’s
ratio $\nu=0.49$. Thus, the shear modulus of the substrate is
$G_{B}={\frac{E_{B}}{2(1+\nu)}}=0.84$ MPa. Scales were modeled to be rigid
with respect to the substrate by imposing rigid body constraints. The contact
between the rigid scales is modeled using the ’Surface to Surface Contact’
algorithm of the software [31]. This contact option is typically used in
literature for biomimetic scale problems of this type [14, 17, 18, 21].

The mechanical loads for bending and twisting were applied quasi-statically to
the system and as boundary conditions. In our numerical studies, we set
$\kappa_{2}=0$. That is, we do not explore bending in the transverse direction
to the substrate, which is considered trivial in the current context. To
compare the combined effect of bending (with strain $\kappa_{1}$) and twisting
(with strain $\kappa_{3}$) with the semi-analytical model presented above, we
perform the simulation in two static steps with nonlinear geometry option
(NLGEOM on) [31]. In the first step, the bending rotations with the same
magnitude but opposite directions were applied to the cross-sections at either
ends. The magnitude of bending was increased linearly from 0 to
(approximately) $\kappa_{1}L_{B}/2$, where $L_{B}$ is the length of the
substrate. The final value could be attained only approximately as edge
effects precluded a constant bending strain could not be maintained throughout
the substrate. To ameliorate these effects we extracted an averaged strain in
the beam far from the ends. More details on this procedure will be explained
in A. In the second step, the bending rotations was fixed at the final value
($\kappa_{1}L_{B}/2$) and the twisting rotations were applied to the both end
cross-sections with reverse directions, again linearly increasing from 0 to
(approximately) $\kappa_{3}L_{B}/2$ during the step time. A frictionless
surface-to-surface contact was applied to the scales surfaces. To test the
reliability of the numerical results, we carried out a mesh convergence study.
Sufficient mesh density was found for different regions of the model
confirming computational accuracy. Mesh convergence led to a total number of
about 230,000 elements. Because of the complex geometry of the system
particularly due scales inclusions in the substrate, the top layer of
substrate was meshed with tetrahedral quadratic elements C3D10 and other
regions were meshed with quadratic hexahedral elements C3D20.

##  4 Results and Discussion

In this section we present results exploring cross-curvature effects, with
various types of scales contact color-coded differently in the plots. Without
loss of generality, we take the initial scale angle $\theta_{0}=0^{\circ}$
(grazing scales) since we are more interested in scale engagement behavior. A
$\theta_{0}>0$ would only shift the plots forward till engagement occurs, with
the curves thereafter following the same trajectories [14, 21]. We use the
following non-dimensional parameters: $\eta=l/d$, $\beta=b/d$, and
$\hat{\kappa}_{i}=\kappa_{i}d$, $i=1,2,3$ to explore the effects.

A typical $\hat{\kappa}_{3}$ (twist) versus $\theta$ (scale angle) behavior
(for a given $\hat{\kappa}_{1}$) is shown in Fig. 6. In this particular
figure, we have used $\hat{\kappa}_{1}=0.23$, $\eta=3$, $\beta=1.25$ and
$\alpha=30^{\circ}$. The presence of an initial bending can lead to additional
scale contact types, not observed for pure twisting investigated previously
[21, 22]. These various contact regimes are color coded in the twist-angle
plot as follows: the black curve is the contact type where the top edge of the
first scale makes contact with the left edge of the second, the red curve
represent configurations where the corner of the first scale makes contact
with the face of the second scale, green curve the configuration where the top
edge of the first scale makes contact with the right edge of the second scale,
and the blue curve represents the configuration where the top edges of the two
scales make contact. We also provide animations of the contact kinematics (see
supplementary material) to aid in better visualization of this system.

Figure 6: $\hat{\kappa}_{3}$ vs $\theta$ with different modes of contact. Here
the fixed parameters are: $\hat{\kappa}_{1}=0.23$, $\eta=3$, $\beta=1.25$ and
$\alpha=30^{\circ}$.

Recall that in the case of pure twisting [21] (or pure bending [14]) the
kinematic locking is typically identified in the strain(curvature/twist)-angle
plot at the point where the slope goes to infinity. However, in the case of
combined load, at locking, the slopes changes sign abruptly from positive to
negative. This critical point is indicated by a blue triangle in Fig. 6 where
the curve is predicted to transition from green to red. However, in reality
configurations beyond this critical point are physically impossible without
violating interpenetration condition. This critical point is thus interpreted
as locking in this case. This non-orthogonal condition of locking is a unique
artifact of combined load system, not found in pure bending or twisting [14,
21]. In the subsequent analysis, we treat post-locking configurations as
forbidden. For the parameter values explored in this work, we find that the
blue curves (which correspond to configurations with top edge contact) always
occur as locked states. Hence, they are not explicitly shown in the figures
below. In addition, contact regimes that correspond to very high twisting
strains are also not considered as they are of little practical significance
and difficult to achieve without introducing significant material
nonlinearity.

We now explore in detail the cross coupling effects of one curvature over the
other. In Fig. 7 (a) we show the effect of bending ($\hat{\kappa}_{1}$) on the
twisting ($\hat{\kappa}_{3}$) versus scale angle ($\theta$). The various
values of $\hat{\kappa}_{1}$ are shown beside their corresponding curves. This
figure is fundamentally different from pure bending [14] or twisting [21]
kinematics. Here, there are distinct regions of kinematics and they are
dependent on the existing bending curvature, with distinct contact regimes
emerging throughout engagement. Thus existence of bending can result in
entirely new type of kinematic behaviors in twisting. The color coded plots
are an indication that there would be kinks in twisting rigidity as contact
regimes between scales undergo change as depicted, for instance in Fig. 6,
inset C to D. This plot can also be used to note the sensitivity of twisting
to existing bending strains. We find that the sensitivity is quite high and
thus any unintended bending either due to sample processing or loading
asymmetries can drastically change the overall kinematics. Another remarkable
feature is that, locking angle for the scale seems to be insensitive to
bending. Hence, the locking envelopes computed from pure twisting would still
hold in this case. Overall, bending strains have a differential effect on the
twisting kinematics, affecting the kinematics trajectory substantially while
leaving the locked state unchanged.

|
---|---
(a) | (b)
|
(c) | (d)
Figure 7: (a) $\theta$ versus $\hat{\kappa}_{3}$ for different values of
bending strains ($\hat{\kappa}_{1}$). (b) $\theta$ versus $\hat{\kappa}_{1}$
for different values of twisting strains ($\hat{\kappa}_{3}$). (c) FE
simulations of $\theta$ versus $\hat{\kappa}_{3}$ for different values of
bending strains ($\hat{\kappa}_{1}$) (d) FE simulations of $\theta$ versus
$\hat{\kappa}_{1}$ for different values of twisting strains
($\hat{\kappa}_{1}$), $\eta=3$, $\beta=1.25$ and $\alpha=30^{\circ}$.

In the same spirit, in Fig. 7 (b), we highlight the effect of twisting on the
bending behavior of the biomimetic scale beam. Here, we find that a twist in
the system also has a significant impact on the bending behavior. With
increasing positive twist, the bending engagement occurs earlier, whereas the
opposite occurs with negative twist. This trend reflects the fact that
twisting on one side has an ‘opening’ effect on the scale, where they move
apart from each other compared to the other side which leads to a ’closing’
effect. The slopes also change significantly indicating a potential effect on
the overall stiffness of the system. The abrupt changes correspond to sudden
changes in contact regimes and shown with different colors. Here again,
interestingly the presence of twist does not affect the locking angle,
analogous to the twisting case.

The validation of these analytical results using FE models are discussed next.
Fig. 7 (c) and (d) are the counter parts of Fig. 7 (a) and (b), respectively.
In Fig. 7 (c) we compare the FE results (shown as dots) with the results of
our model (shown as dashed lines) for various values of pre-set bending
strains ($\hat{\kappa}_{1}$) in the substrate. The kinks seen in the curves
(close to $\hat{\kappa}_{3}=0.05$) are precisely the transition from red to
green curves shown in Fig. 7 when the mode of contact changes. Analogous
results for $\theta$ versus $\hat{\kappa}_{1}$ for different values of pre-set
twisting strains $\hat{\kappa}_{3}$ are shown in Fig. 7 (d). Note that the
pre-set strains for the two set of simulations approximate. This is because
when the second quasi-static step (see Sec. 3) is performed, the strain values
fixed in during the first quasi-static step does not remain fixed. That is, as
we twist a pre-bent beam, the bending curvatures change mildly under the
twisting deformation. In Fig. 7 (c) and (d), we account for this change by
interpolating over the changing values of $\hat{\kappa}_{1}$ and
$\hat{\kappa}_{3}$, respectively.

We now use our model to explore the twisting kinematics at a given bending
strain, with changing geometric parameters, $\alpha$ and $\eta$. These
parameters have shown to be critical in dictating the overall kinematics of
pure twisting [21]. For this, we fix the bending curvature at
$\hat{\kappa}_{1}=0.1$. First, we probe the significance of overlap ratio
$\eta$ on the overall kinematics of the system. In Fig. 8 (a), we plot
$\theta$ vs $\hat{\kappa}_{3}$ for different values of $\eta$ (for
$\beta=1.25$, $\alpha=30^{\circ}$ and $\kappa_{1}=0.15$). Here, we first note
the similarities with pure twisting case. Like pure twisting, increasing
$\eta$ leads to steeper slopes in the $\theta-\hat{\kappa}_{3}$ plots.
However, the twist completely changes the kinematic trajectories by
introducing two distinct contact regimes. The abrupt changes in contact are
more apparent for lower overlap ratio but are ameliorated at higher overlaps.
Note that on the right side, the plots continue to increase in steepness until
the locked states are reached. On the other hand, on the left side (negative
twist), the decrease in slope (flatten) till they collapse into a point. This
is the point of disengagement ($\kappa_{3}=\kappa_{3}^{*}$) of the scales
beyond which, they no longer remain in contact. A surprising outcome of this
plot is that, the disengagement point seems to be independent of the overlap
ratio. This can be understood if this point is visualized for an RVE, Fig. 9
(intersection between the two scales is shown in red). Evidently, at this
point, the nature of contact is the entire edge of one scale moving over the
adjacent one. Thus, at $\kappa_{3}^{*}$, $\eta$ itself is not well defined and
the scale angle is purely determined by the dihedral angle of two adjacent
intersecting planes that is dependent only on the applied strains.

We see here again that $\eta$ influences the angle of locking. We find that
this parameter has a significant role to play in determining the overall
kinematics of the system. Increasing $\eta$ leads to an overall stiffer
nonlinearity but at the same time can potentially change the overall contact
regime sequence up until locking.

|
---|---
(a) | (b)
|
(c) | (d)
Figure 8: (a) $\theta$ versus $\hat{\kappa}_{3}$ for different values of
$\eta$ (for $\beta=1.25,\;\alpha=30^{\circ},\;\hat{\kappa}_{1}=0.1$). (b)
$\theta$ versus $\kappa_{3}$ for different values of $\alpha$
($\eta=3,\;\beta=1.25$, $\kappa_{1}=0.1$). (c) $\theta$ versus $\kappa_{1}$
for different values of $\eta$ (for
$\beta=1.25,\;\alpha=30^{\circ},\;\kappa_{3}=0.1$). (d) $\theta$ versus
$\kappa_{1}$ for different values of $\alpha$ ($\eta=3,\;\beta=1.25$,
$\kappa_{3}=0.1$).

Next, In Fig. 8 (b), we plot $\theta$ vs $\hat{\kappa}_{3}$ for different
values of $\alpha$ (for $\beta=1.25$, $\eta=3$ and $\hat{\kappa}_{1}=0.1$).
The effect of the tilt angle $\alpha$ is dramatic. For relatively small
angles, i.e. $\alpha<20^{\circ}$, there is a relatively ‘stiff’ response on
either direction of twisting. These reflect scales sliding on either direction
of twisting. However, there is a marked anisotropy between the directions in
both magnitude and contact regimes. More interestingly, at higher tilt angles
such bi-directionally distinct stiffness disappears altogether indicating the
scales do not or only weakly engage in other other direction. The lack of
engagement in the other direction can be visualized as ‘opening’ of the scales
in one direction vs closing. The effect of $\alpha$ on locking is also
pronounced. At lower tilt angles, the stiffening effect disappears into a
lockless behavior. This is due to scales sliding past each other without ever
satisfying the contact constraints. Locking behavior emerges again for higher
tilt angles.

Figure 9: Visual depiction of the scales in contact at the point of
disengagement as scales ’open’ up at $\kappa_{3}=\kappa_{3}^{*}$. (for
$\beta=1.25,\;\alpha=30^{\circ},\;\hat{\kappa}_{1}=0.1$)

In Fig. 8 (c), we now look at the sensitivity of the bending kinematics in the
presence of twist. For these plots the we keep the geometry of the substrates
the same and now impose a twist of $\hat{\kappa}_{3}=0.1$ before bending and
keep the tilt angle $\alpha=30^{\circ}$. Here, we see that the overall impact
of higher $\eta$ is to increase the slope of $\theta-\hat{\kappa}_{1}$ plot,
similar to pure bending [14]. However, the presence of twist changes the
nature of the $\theta-\hat{\kappa}_{1}$ curve by introducing different contact
regimes, very similar to the twisting case discussed earlier (Fig. 8 (a)).
Similar to that case, the abruptness of the kinks diminish at higher $\eta$,
disappearing altogether when $\eta=5$. Thus, higher scale overlap has a
suppressing effect on contact regime transitions. Interestingly, the scales
remain engaged even at a negative bending value, due to the ‘closing’ effect
on scales due to twisting where scales are pushed closer to each other.
However, eventually the scales again lose contact as bending on the other
direction increases sufficiently, once again exhibiting the ’convergence’ at
the point of disengagement described previously.

We also investigate the effect of an existing twist on
$\theta-\hat{\kappa}_{1}$ relationship for various tilt angles $\alpha$ but
with fixed $\eta=3$ in Fig. 8 (d). Here, we find that at lower tilt angle, the
contact regime transition is absent, appearing as tilt angle increases. This
shows that the tilt angle is an important parameter that controls the
complexity of contact regimes. Overall, the effect of pre-twist in bending is
muted at lower tilt angles and higher overlap and more pronounced otherwise.

The investigation above, set the stage to investigate the positive-negative
effects of twist-bend combinations. To put this bend-twist kinematics in more
general context, we compute the curvature dependence of scale angles for a
biomimetic scale beam of fixed geometric parameters. This results in a phase
plot as shown in Fig. 10. The curvature limits are taken to maintain self
consistency with respect to contact regimes studied. For higher twist levels,
more contact possibilities would need to be included. The phase plot in Fig.
10 is by solving the kinematic relationships developed earlier connecting
scale angles and curvatures numerically at various values of bending and
twisting strains.

Figure 10: Slice of the engagement region, $\Gamma_{e}$, as defined in (21)
for $\hat{\kappa}_{2}=0$. Points lying inside the shaded region correspond to
configurations where scales are in contact. The color bar shows the scale
angle $\theta$ (in radians) for the corresponding configuration. Here
$\theta_{0}=5^{\circ}$ and $\alpha=30^{\circ}$, $\eta=3$, and $\beta=1.25$.

The phase boundary on the left is the limit of engagement, i.e. we do not have
any scales engagement beyond this boundary. The phase boundary on the right
hand side is the locking boundary. Any $(\hat{\kappa}_{1},\hat{\kappa}_{3})$
combination is forbidden in this region. Thus, the region in between indicates
the region of engagement. It is interesting to note that there is no inherent
symmetry between bending and twisting. In other words, the effect of bending
on twisting is fundamentally different from the converse. In addition, one can
see that multiple bend-twist combinations can give rise to same scale angles.
This angular degeneracy is quite remarkable and could have significance for
inverse designs.

Finally, we discuss the mechanical behavior of these systems. In Fig. 11 (a)
we plot bending moment vs bending strains and discover the role played by the
presence of twist for a given $\eta$. Higher values of positive twist shift
the engagement to earlier parts, whereas higher values of negative twist shift
it in the opposite direction. Interestingly, unlike smooth plots seen for pure
bending earlier [14], twist effectively changes the contact regime in terms of
discontinuities and jumps in the plot indicating sudden changes in bending
rigidity. When a twist is added in the positive direction, the neutral
position (no bending strain) is fully engaged (due to the engagement brought
about by the twist). The FE comparison with these plots is shown in Fig. 11
(b). Notice the validation region is for much smaller values of bending and
twisting. This is due to inherent limitation of traditional FE software in
simulating this type of system and already well known in literature [14, 21,
32]. The agreements are excellent with respect to the theoretical values. The
observed deviations are expected both due to the complex nature of contact,
imposing global periodicity, and also the difficulties in keeping one
curvature constant in an actual numerical simulation due to edge effects.

|
---|---
(a) | (b)
Figure 11: Bending moment $M_{1}$ versus $\hat{\kappa}_{1}$ for various values
of $\hat{\kappa}_{3}$: (a) Plots for the analytical model (b) Plots using FE
analysis with dashed lines showing theoretical predictions.

We see similar and even more dramatic effect in the twist-torque diagram in
Fig. 12 (a). Here again, we see the role of the bending strains as the torque-
twists which are smooth nonlinear plots in pure twisting case [21] now turn
into highly discontinuous plots with disparate twist modulus. Yet again, an
existing bending strain can cause the neutral position (zero twist strain) to
be pre-engaged due to bending. These highlight the appreciable differences
that can be brought about from cross-coupling effects. Fig. 12 (b) shows the
companion FE torque-twist plots that also indicate excellent agreement with
the analytical results. As in the case for bending, the deviations are mainly
because of the difficulties in keeping the cross curvatures fixed for the
duration of the simulation. The other sources of errors are possibly due to
assuming fixed curvatures as opposed to spatially varying curvatures arising
in a FE simulation. Analogous to bending case, imposing global periodicity in
scale engagements also leads to some differences in behavior. Please note that
many of these non ideal effects and sources of discrepancies have been
discussed earlier [16, 17].

|
---|---
(a) | (b)
Figure 12: Twisting moment $M_{3}$ versus $\hat{\kappa}_{3}$ for various
values of $\hat{\kappa}_{1}$: (a) Plots for the analytical model (b) Plots
using FE analysis with dashed lines showing theoretical predictions.

##  5 Conclusions

In this work, we addressed the cross coupling effects of bending and twisting
in a biomimetic scale elastic beam for the first time. Here the scales were
plate like rectangular inclusions protruding at an angle from the surface of
the elastic substrate. We find highly intricate and often surprising effect of
one over the other across the kinematics and mechanics. We quantified these
effects by developing analytical relationships within the framework of
Cosserat kinematics and global-local energy balance. This model reduced to the
earlier developed model for pure bending and twisting in literature and was
also validated with FE simulations. This study completes a significant missing
piece in the mechanics of fish scale inspired biomimetic system, which is of
great practical importance in applications ranging from high performance smart
skins to soft robotic systems.

## Acknowledgement

This work was supported by the United States National Science Foundation’s
Civil, Mechanical, and Manufacturing Innovation, CAREER Award #1943886.

## References

  * [1] B. Wang, W. Yang, V. R. Sherman, M. A. Meyers, Pangolin armor: overlapping, structure, and mechanical properties of the keratinous scales, Acta biomaterialia 41 (2016) 60–74.
  * [2] I. H. Chen, J. H. Kiang, V. Correa, M. I. Lopez, P.-Y. Chen, J. McKittrick, M. A. Meyers, Armadillo armor: mechanical testing and micro-structural evaluation, Journal of the mechanical behavior of biomedical materials 4 (5) (2011) 713–722.
  * [3] C. LaTorre, B. Bhushan, Investigation of scale effects and directionality dependence on friction and adhesion of human hair using afm and macroscale friction test apparatus, Ultramicroscopy 106 (8-9) (2006) 720–734.
  * [4] K. Michielsen, D. Stavenga, Gyroid cuticular structures in butterfly wing scales: biological photonic crystals, Journal of The Royal Society Interface 5 (18) (2008) 85–94.
  * [5] M. J. Buehler, Nature designs tough collagen: explaining the nanostructure of collagen fibrils, Proceedings of the National Academy of Sciences 103 (33) (2006) 12285–12290.
  * [6] U. G. Wegst, H. Bai, E. Saiz, A. P. Tomsia, R. O. Ritchie, Bioinspired structural materials, Nature materials 14 (1) (2015) 23–36.
  * [7] K. Kertész, G. Molnár, Z. Vértesy, A. Koós, Z. Horváth, G. Márk, L. Tapasztó, Z. Bálint, I. Tamáska, O. Deparis, et al., Photonic band gap materials in butterfly scales: A possible source of “blueprints”, Materials Science and Engineering: B 149 (3) (2008) 259–265.
  * [8] J. Long, M. Hale, M. Mchenry, M. Westneat, Functions of fish skin: flexural stiffness and steady swimming of longnose gar, lepisosteus osseus, Journal of Experimental Biology 199 (10) (1996) 2139–2151.
  * [9] J. Song, S. Reichert, I. Kallai, D. Gazit, M. Wund, M. C. Boyce, C. Ortiz, Quantitative microstructural studies of the armor of the marine threespine stickleback (gasterosteus aculeatus), Journal of structural biology 171 (3) (2010) 318–331.
  * [10] S. H. Sadati, Y. Noh, S. E. Naghibi, A. Kaspar, T. Nanayakkara, Stiffness control of soft robotic manipulator for minimally invasive surgery (mis) using scale jamming, in: International Conference on Intelligent Robotics and Applications, Springer, 2015, pp. 141–151.
  * [11] Y. Wei, Y. Chen, T. Ren, Q. Chen, C. Yan, Y. Yang, Y. Li, A novel, variable stiffness robotic gripper based on integrated soft actuating and particle jamming, Soft Robotics 3 (3) (2016) 134–143.
  * [12] E. T. Roche, M. A. Horvath, I. Wamala, A. Alazmani, S.-E. Song, W. Whyte, Z. Machaidze, C. J. Payne, J. C. Weaver, G. Fishbein, et al., Soft robotic sleeve supports heart function, Science Translational Medicine 9 (373) (2017) eaaf3925.
  * [13] J.-Y. Sire, P. C. Donoghue, M. K. Vickaryous, Origin and evolution of the integumentary skeleton in non-tetrapod vertebrates, Journal of Anatomy 214 (4) (2009) 409–440.
  * [14] R. Ghosh, H. Ebrahimi, A. Vaziri, Contact kinematics of biomimetic scales, Applied Physics Letters 105 (23) (2014) 233701.
  * [15] R. Ghosh, H. Ebrahimi, A. Vaziri, Frictional effects in biomimetic scales engagement, EPL (Europhysics Letters) 113 (3) (2016) 34003.
  * [16] R. Ghosh, H. Ebrahimi, A. Vaziri, Non-ideal effects in bending response of soft substrates covered with biomimetic scales, Journal of the mechanical behavior of biomedical materials 72 (2017) 1–5.
  * [17] H. Ali, H. Ebrahimi, R. Ghosh, Bending of biomimetic scale covered beams under discrete non-periodic engagement, International Journal of Solids and Structures 166 (2019) 22–31.
  * [18] H. Ali, H. Ebrahimi, R. Ghosh, Tailorable elasticity of cantilever using spatio-angular functionally graded biomimetic scales, Mechanics of Soft Materials 1 (2019) 10.
  * [19] H. Ali, H. Ebrahimi, R. Ghosh, Frictional damping from biomimetic scales, Scientific Reports 9 (1) (2019) 14628.
  * [20] H. Ali, H. Ebrahimi, J. Stephen, P. Warren, R. Ghosh, Tailorable stiffness lightweight soft robotic materials with architectured exoskeleton, in: AIAA Scitech 2020 Forum, 2020, p. 1551.
  * [21] H. Ebrahimi, H. Ali, R. A. Horton, J. Galvez, A. P. Gordon, R. Ghosh, Tailorable twisting of biomimetic scale-covered substrate, EPL (Europhysics Letters) 127 (2) (2019) 24002.
  * [22] H. Ebrahimi, H. Ali, R. Ghosh, Coulomb friction in twisting of biomimetic scale-covered substrate, Bioinspiration & Biomimetics 15 (5) (2020) 056013.
  * [23] H. Ebrahimi, H. Ali, J. Stephen, S. Dharmavaram, R. Ghosh, Emergent mechanical properties of biomimetic exoskeletal metamaterials, in: Bioinspiration, Biomimetics, and Bioreplication XI, Vol. 11586, International Society for Optics and Photonics, 2021, p. 115860N.
  * [24] H. Ebrahimi, H. Ali, J. Stephen, R. Ghosh, Fish scales: Primitive basis for modern metamaterials, EPL (Europhysics Letters) 133 (6) (2021) 68001.
  * [25] F. J. Vernerey, F. Barthelat, On the mechanics of fishscale structures, International Journal of Solids and Structures 47 (17) (2010) 2268–2275.
  * [26] F. J. Vernerey, F. Barthelat, Skin and scales of teleost fish: Simple structure but high performance and multiple functions, Journal of the Mechanics and Physics of Solids 68 (2014) 66–76.
  * [27] J. Kemoole, (accessed September 5, 2020) (2017).  [link].
URL https://flic.kr/p/YSe1Us

  * [28] S. S. Antman, Nonlinear Problems of Elasticity, Springer, 1995.
  * [29] B. Audoly, Y. Pomeau, Elasticity and geometry: from hair curls to the non-linear response of shells, Oxford university press, 2010.
  * [30] A. C. Ugural, S. K. Fenster, Advanced mechanics of materials and applied elasticity, Pearson Education, 2011.
  * [31] V. Abaqus, 6.14 documentation, Dassault Systemes Simulia Corporation 651 (6.2) (2014).
  * [32] A. Shafiei, J. W. Pro, R. Martini, F. Barthelat, The very hard and the very soft: Modeling bio-inspired scaled skins using the discrete element method, Journal of the Mechanics and Physics of Solids 146 (2021) 104176.

##  Appendix A Extracting Cosserat Strains from FE Data

To compare the FE simulations with the analytical model presented in Sec. 2,
we extract the Cosserat bending strains $\kappa_{1}$, $\kappa_{2}$, and
$\kappa_{3}$ from the FE data through the following steps. First, we extract
the position vector for the mid-point of the top face of the beam from the
numerical simulations. This vector is an estimate for $\mathbf{r}(s)$ as given
in (8). The information of the directors $\mathbf{d}_{1}(s)$ and
$\mathbf{d}_{2}(s)$ is extracted by subtracting the position vectors of
center-line and the right-edge of the beam, respectively, from the estimate
for $\mathbf{r}(s)$ and normalizing the resulting quantities to produce unit
vectors. After validating inextensibility by verifying that the change in
length of the rod is small (with typical percentage relative error $\Delta
L_{B}/L_{B}\approx 1\%$), we compute the third director using
$\mathbf{d}_{3}=\mathbf{r}^{\prime}(s)$, estimating the right hand side using
finite differences along the beam. It follows from (1) that rotation matrix
mapping the cross-sections of the beam is given by
$\mathbf{Q}(s)=[\mathbf{d}_{1}(s),\mathbf{d}_{2}(s),\mathbf{d}_{3}(s)]$ where
the directors are taken to be the column vectors of the matrix. The skew-
symmetric matrix $\mathbf{K}$ containing the bending strains
$\kappa_{1}(s),\;\kappa_{2}(s),\;\kappa_{3}(s)$ (cf. (4)) along the length of
the beam is computed using the formula
$\mathbf{K}=\mathbf{Q}^{\prime}(s)\mathbf{Q}^{-1}(s)$, where
$\mathbf{Q}^{\prime}(s)$ is estimated using finite differences. To avoid
boundary effects, we average the bending strains over the middle half of the
beam to obtain estimates for average bending strains in the beam.

Figure 13: Extracting Cosserat strain from FE Data.


# The functional equation for $\zeta$

Keith Ball

This note contains a short proof of the functional equation for the zeta
function. The argument is computationally similar to the Hankel contour
arguments found in many places but is made simpler to visualise by the
symmetries of the sine. In this form the argument seems to be new or at least
not widely known. We begin with a simple lemma which uses integration by parts
to introduce the sine (or hyperbolic sine).

######  Lemma 1.

For $s>1$

$\Gamma(s)\zeta(s)=\frac{1}{4s}\int_{0}^{\infty}\frac{t^{s}}{\sinh^{2}(t/2)}\,dt.$

Proof

$\displaystyle\Gamma(s)\zeta(s)$ $\displaystyle=$
$\displaystyle\sum_{1}^{\infty}\frac{1}{n^{s}}\int_{0}^{\infty}e^{-t}t^{s-1}\,dt=\sum_{1}^{\infty}\int_{0}^{\infty}e^{-nt}t^{s-1}\,dt$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\frac{t^{s-1}}{e^{t}-1}\,dt=\frac{1}{s}\int_{0}^{\infty}\frac{e^{t}\,t^{s}}{(e^{t}-1)^{2}}\,dt$
$\displaystyle=$
$\displaystyle\frac{1}{4s}\int_{0}^{\infty}\frac{t^{s}}{\sinh^{2}(t/2)}\,dt.$

We now move on to the theorem itself.

######  Theorem 2.

The function $s\mapsto(s-1)\zeta(s)$ defined on $(1,\infty)$ has an analytic
continuation to the entire complex plane and for all $s$ we have

$\zeta(s)=2(2\pi)^{s-1}\sin(\pi s/2)\Gamma(1-s)\zeta(1-s)$

with the usual convention regarding the removable singularities.

Proof Consider the integral

$\frac{1}{2\pi i}\int_{C}\frac{\pi^{2}z^{1-s}}{\sin^{2}\pi z}\,dz$

where the contour $C$ is a vertical line between 0 and 1, say $(\Re z=1/2)$,
traversed upwards. In this and the rest of the argument we take powers
$z^{-s}$ only for complex numbers satisfying $\Re z\geq 0$ so there is no
ambiguity. Because of the rapid growth of $\sin\pi z$ as $z$ moves away from
the real axis, this function is easily seen to be an entire function of $s$.

If $s>1$ then the factor $(x+iy)^{1-s}\rightarrow 0$ uniformly in $y$ as
$x\rightarrow\infty$ while the factor $\sin^{2}\pi(x+iy)$ is periodic in $x$,
so we can compute the integral as the negative of the sum of the residues of

$\frac{\pi^{2}z^{1-s}}{\sin^{2}\pi z}$

at the positive integers. The residue at $n$ is $(1-s)/n^{s}$ so we find that
the integral is $(s-1)\zeta(s)$.

On the other hand, if $s<0$ the integral makes sense at $z=0$ so we may shift
the contour onto the imaginary axis. By splitting into the upper and lower
halves of the axis we get that the integral is

$\displaystyle\frac{(e^{i\pi/2})^{2-s}-(e^{-i\pi/2})^{2-s}}{2\pi
i}\int_{0}^{\infty}\frac{\pi^{2}y^{1-s}}{\sin^{2}\pi iy}\,dy$ $\displaystyle=$
$\displaystyle-\pi\sin(\pi s/2)\int_{0}^{\infty}\frac{y^{1-s}}{\sinh^{2}\pi
y}\,dy$ $\displaystyle=$ $\displaystyle\frac{-\pi\sin(\pi
s/2)}{(2\pi)^{2-s}}\int_{0}^{\infty}\frac{t^{1-s}}{\sinh^{2}(t/2)}\,dt$
$\displaystyle=$ $\displaystyle\frac{-\pi\sin(\pi
s/2)}{(2\pi)^{2-s}}4(1-s)\Gamma(1-s)\zeta(1-s)$ $\displaystyle=$
$\displaystyle 2(2\pi)^{s-1}\sin(\pi s/2)(s-1)\Gamma(1-s)\zeta(1-s)$

where the last but one identity uses the lemma.

Keith Ball
Mathematics Institute
University of Warwick
Coventry CV4 7AL

<EMAIL_ADDRESS>