We present a model system for strongly nonlinear transition waves generated in a periodic lattice of bistable members connected by magnetic links. The asymmetry of the on-site energy wells created by the bistable members produces a mechanical diode that supports only unidirectional transition wave propagation with constant wave velocity. We theoretically justify the cause of the unidirectionality of the transition wave and confirm these predictions by experiments and simulations. We further identify how the wave velocity and profile are uniquely linked to the double-well energy landscape, which serves as a blueprint for transition wave control

Arrieta, Andres F.

Chong, Christopher

Daraio, Chiara

Kochmann, Dennis M.

Nadkarni, Neel

Caltech Authors - Main

SUPPLEMENTAL MATERIAL for:
Unidirectional Transition Waves in Bistable Lattices
by N. Nadkarni, A. F. Arrieta, C. Chong, D. M. Kochmann, C. Daraio
1 Snapshot of the transition wave
Figure 1: Snap-shot sequence showing the transition wave as it propagates through the experimental
lattice. Images were acquired at 4000 fps. As can be seen, the wave is highly localized to a few
elements.
2 Theoretical derivation of the continuum model
Consider the infinite nonlinear lattice with on-site (bistable) and inter-site nonlinearities with gov-
erning equation,
mun,tt +A(un+1 − un + L)p −A(un − un−1 + L)p + αun,t + βφ′(un) = 0. (1)
where L is the lattice parameter and p < −1. We can therefore deduce the following relations:
1. The mass density should be a constant: m/L = constant ⇒ m = O(L).
2. The energy density of the nonlinear spring should be constant: − Ap+1 (un+1 − un + L)p+1 /L =
− Ap+1 (L+ L)p+1 /L = constant ⇒ A = O(L−p) where  ∼ O(1).
3. The energy density of the bistable function should be constant: βφ′(un)/L = constant ⇒ β =
O(L).
4. The dissipation potential density should be a constant: 12αu
2
n,t/L = constant ⇒ α = O(L).
Assuming a traveling wave solution of the form, un(t) = u(nL− vt) = u(ξ) gives
mv2uξξ +A (u(ξ + L)− u(ξ) + L)p −A (u(ξ)− u(ξ − L) + L)p − vαuξ + βφ′(u) = 0. (2)
Using appropriate Taylor expansions for u(ξ + L) and u(ξ − L) results in
mv2uξξ +AL
p
(
1 + uξ +
L
2 uξξ + ...
)p −ALp (1 + uξ − L2 uξξ + ...)p − vαuξ + βφ′(u) = 0. (3)
If we analyze the inter-element forcing function in the continuum limit, we obtain
F = −ALp (1 + uξ + L2 uξξ + ...)p . (4)
1
The negative sign is required because the external force is compressive. As L→ 0, ALp → k ∼ O(1),
giving,
F = −k (1 + uξ)p (5)
As p < −1, for (5) to make physical sense, we need to the force to be finite which implies that
1 + uξ 6= 0. Referring back to the governing equation, we can perform another Taylor expansion in
the continuum limit (L→ 0) so that
mv2uξξ +AL
p
(
(1 + uξ)
p + L2 p uξξ(1 + uξ)
p−1 + ...
)−ALp ((1 + uξ)p − L2 p uξξ(1 + uξ)p−1 + ...)
− vαuξ + βφ′(u) = 0.
(6)
Dividing throughout by L and rearranging gives
ρv2uξξ +
ALp+1
L
(
p uξξ(1 + uξ)
p−1)+ ALp+3L [16(p− 1)puξξuξξξ(uξ + 1)p−2+
1
24(p− 2)(p− 1)pu3ξξ(uξ + 1)p−3 + 112puξξξξ(uξ + 1)p−1]− v αLuξ + βLφ′(u) = 0.
(7)
We can see that, m/L = ρ ∼ O(1), −ApLp = ρc20 ∼ O(1), α/L = γ ∼ O(1) and β/L ∼ O(1).
Define βLφ
′(u) = ψ′(u). Taking the limit with respect to L (while keeping O(L2) terms), we obtain
the continuum equation
ρv2uξξ − ρc20 uξξ(1 + uξ)p−1 − 124L2ρc20[4(p− 1)uξξuξξξ(uξ + 1)p−2
+ (p− 2)(p− 1)u3ξξ(uξ + 1)p−3 + 2uξξξξ(uξ + 1)p−1]− vγuξ + ψ′(u) = 0.
(8)
Multiplying by uξ and integrating from −∞ to ∞, we obtain
− 124L2ρc20
∫ ∞
−∞
[4(p− 1)uξξuξξξ(uξ + 1)p−2 + (p− 2)(p− 1)u3ξξ(uξ + 1)p−3 + 2uξξξξ(uξ + 1)p−1]uξdξ
+
∫ ∞
−∞
(
ρv2 − ρc20 (1 + uξ)p−1
)
uξuξξdξ +
∫ ∞
−∞
ψ′(u)uξdξ = vγ
∫ ∞
−∞
u2ξ dξ.
(9)
Now, as the system is dissipative, it is fair to assume that, as t→∞ or ξ → −∞, all derivatives of
u go to zero. Also, as the system is initially at rest, all derivatives vanish as t → −∞ or ξ → ∞.
With this in mind, we review the second integral in (9). If the integration variable is rewritten
using uξξdξ = d(uξ) and w = uξ, we get∫ ∞
−∞
(
ρv2 − ρc20 (1 + uξ)p−1
)
uξuξξdξ =
∫ w(ξ→∞)
w(ξ→−∞)
(
ρv2 − ρc20 (1 + w)p−1
)
wdw = 0, (10)
as 1 + w 6= 0, as proved earlier. Next, let us analyze the second term in the first integral (without
the constant − 124L2c20) in (9) and integrate by parts,∫ ∞
−∞
[(p− 1)(p− 2)(uξ + 1)p−3uξξ][u2ξξuξ]dξ =− 2
∫ ∞
−∞
(p− 1)(uξ + 1)p−2uξξuξuξξξdξ
−
∫ ∞
−∞
(p− 1)(uξ + 1)p−2u3ξξdξ.
(11)
Performing a similar integration by parts on the last term of the first integral in (9), we obtain
2
∫ ∞
−∞
uξξξξ(uξ + 1)
p−1uξdξ =− 2
∫ ∞
−∞
(p− 1)(uξ + 1)p−2uξuξξuξξξdξ
+
∫ ∞
−∞
(p− 1)(uξ + 1)p−2u3ξξdξ.
(12)
2
Therefore, all the terms of the first integral will cancel each other. Similarly, if u(ξ →∞) = ui and
u(ξ → −∞) = uf , and without loss of generality, if v > 0, then (9) becomes∫ ∞
−∞
ψ′(u)uξdξ =
∫ ψ(ui)
ψ(uf )
d(ψ(u)) = ψ(ui)− ψ(uf ) = vγ
∫ ∞
−∞
u2ξ dξ ≥ 0⇒ ψ(ui) ≥ ψ(uf ) (13)
3 Experimental force characterization
The magnetic force and bistable snapping force measurements were made using a Zwick-Roell
compression testing machine. For measuring the magnetic interaction force, one of the magnets
was kept fixed at the base of the machine while the other was fixed to the moving tip of the
machine. A schematic of the experimental setup is shown in Fig. 1. The tip was moved downward
at a constant rate of 20 mm/min and the force on the tip was recorded by the machine. A similar
Magnets 
v 
F 
Computer recording force 
‘F’ and displacement ‘u’ 
u 
Compression testing machine 
Figure 2: Schematic of maximum bistable force measurement
test was conducted for measuring the snapping force of the bistable composite. A schematic of
the experiment is shown in Fig. 2. The element was held fixed, while the magnet on the machine
tip was used to provide a non-contact force so as to allow the element to deform freely. In this
experiment as well, the tip was moved downward at a constant rate of 20 mm/min.
4 Experimental details
A novel composite design featuring a spatially varying fiber distribution is utilized to produce
bistable members with a tailored strain potential topology. This is comprised of a central [0/90]
unsymmetrically laminated section surrounded by two transition sections, as shown in Fig. 4. The
proposed lamination distribution allows for tailoring the depths of the potential wells and associated
snapping forces, while enabling clamping the short edges. The dimensions for the elements used
in the model demonstrator used in the main manuscript are given in table 1. Table 2 provides
the material properties of the used c-m-p (CM-Preg T-C-120/625 CP002 35) prepreg system. It
3
v Computer recording force 
‘F’ and displacement ‘u’
Compression testing machine
Magnets F
Clamps
u
Bistable 
element
Figure 3: Schematic of maximum bistable force measurement
is worth mentioning that the dimensions of the regions making up the fibre distribution of this
laminates can provide a very broad range of different potential wells and snapping forces as detailed
in Ref. [1]. The magnets were of the R-19-09-06 N type (mass of 10 g, inner diameter 9.5 mm,
x
y
L
b
z
x
[0o/0o] [90o/0o] [0o/90o] [90o/0o] [0o/0o]
0o
90o
Runsym,centralBCs
[0o/0o] [0o/0o]
L1 L2 L6 L7L5L4L3
h
BCs
Figure 4: Tailored distribution of used bistable composite laminates.
outer diameter 19.1 mm and thickness of 6.4 mm), supplied by Supermagnete. The experiments
were performed with two Photron Ux 100 cameras, using a Digital Image Correlation system from
Correlated Solutions.
4
L b h L1 L2 L3 L4 L5 L6 L7
[mm] [mm] [mm] [mm] [mm] [mm] [mm] [mm] [mm] [mm]
220 64 0.25 31 15 5.0 120 5.0 15 31
Table 1: Geometric properties of the sections composing the spatially varying fiber for the bistable
elements. Refer to for the schematic representation in Fig. 4 of the given parameters.
Material Fibre vol. E11 E22 G12 ν12 ρ α11 α22
[%] [GPa] [GPa] [GPa] [-] [ kg
m3
] K−1
CFRP 60 161 10 4.4 0.3 1570 -1.8E-8 2.25E-5
Table 2: Material properties for a typical ply of CFRP c-m-p (CM-Preg T-C-120/625 CP002 35)
prepreg used to manufacture the bistable elementsly. Nominal prepreg thickness 0.125 mm.
5 Numerical simulations for exact traveling waves
Transition waves of Eq. (1) (of the manuscript) correspond to orbits of Eq. (2) (of the manuscript)
that connect the equilibrium u = ui and u = uf where ui and uf are the positions of the potential
wells of βφ. Such special orbits of advance-delay equations, like Eq. (2) (of the manuscript), can
be approximated numerically up to a prescribed tolerance [2]. In particular, one can make the
discretization uj := u(j∆ξ) where j ∈ Z which, upon a choice of finite difference approximation of
the derivatives, results in an algebraic equation that can be solved via Newton iterations:
v2m
uj+1 − 2uj + uj−1
∆ξ2
− vαu
j+1 − uj−1
2∆ξ
+ βφ′(uj)
+A(uj+q − uj + L)p −A(uj − uj−q + L)p = 0,
(14)
where ∆ξ is chosen such that q = L/∆ξ is an integer. We impose no flux boundary conditions and
let the wave velocity v be a variable of the system (rather than a fixed parameter). We found that
despite varying the initial guess, our algorithm converges to the same profile and wave velocity v,
implying that for a fixed set of system parameters, there is a unique wave velocity of the transition
wave.
6 Robustness analysis
6.1 Sensitivity of wave velocity with respect to interaction coefficient
Fig. 5 shows the variation of the wave velocity for various coefficients of the interaction potential
(parameter p of Eq. (1) in the manuscript). It can be seen that the velocity becomes zero beyond
a certain critical value of the parameter p. From a physical point of view this is reasonable.The
strength of inter-element forcing decreases with decreasing nonlinearity parameter p. When the
initial element is snapped over, the small force due to magnetic coupling not sufficient for the
subsequent element to snap over the on-site potential energy barrier, thereby causing the wave
to stagnate. This is analogous to the propagation failure observed in discrete reaction-diffusion
5
lattices [3]. The value of p = −3.274 in the experiments lies well above the critical range for each of
the parameter values shown in Fig. 5. For a fixed lattice distance, the potential barrier decreases as
the rail distance increases. Therefore, it is not surprising that the critical value of the nonlinearity
p at which the wave stoppage occurs decreases as the rail distance increases, see Fig. 5(a). For a
constant rail distance, as the lattice distance increases, the interaction force between magnets in
the undeformed state decreases (since the distance between the magnets that couple the elements is
increasing). Hence, the nonlinearity p at which the wave stops propagating increases as the lattice
distance increases, see Fig. 5(b).
interaction potential coefficient (p)
-6 -5 -4 -3 -2 -1
sp
ee
d 
in
 m
/s
0
2
4
6
8
10
R = 21.5cm, L = 6cm
R = 22cm, L = 6cm
R = 22.5cm, L = 6cm
(a)
interaction potential coefficient (p)
-6 -5 -4 -3 -2 -1
sp
ee
d 
in
 m
/s
0
1
2
3
4
5
6
7
8
9
R = 22.5cm, L = 6cm
R = 22.5cm, L = 7cm
R = 22.5cm, L = 8cm
(b)
Figure 5: (a) Variation of wave velocity with interaction coefficient p for different rail distances while
keeping the lattice distance constant. (b) Variation of wave velocity with interaction coefficient p
for different lattice distances while keeping the rail distance constant.
6.2 Sensitivity of wave velocity with respect to asymmetry of the on-site bistable
potential
To obtain a theoretical approximation of the onsite potential ψ(u) for the snapping elements, the
experimentally measured values of the distance and corresponding force were fit by a polynomial
spline (see Fig. 1a of the manuscript). To investigate the role of the asymmetry numerically we
introduce an asymmetry parameter  where  = 0 corresponds to the asymmetry of the experimen-
tally measured on-site potential, and  = 1 corresponds to a symmetric potential (see Fig. 6(a) and
text below for details). As an example, we consider a rail distance of R = 22.5 cm and a lattice
spacing of L = 8 cm. We performed a parametric continuation of roots of Equation (3) of the main
manuscript with respect to the asymmetry parameter . The corresponding value of the wave speed
v is plotted against the asymmetry parameter  in Fig. 6(b). Note the wave speed decreases as the
asymmetry weakens. It appears that there is a threshold wave speed, since solutions cease to exist
for a finite, non-zero value of the asymmetry parameter  (about cr = 0.754). This is also in line
with earlier observations that wave propagation is possible in the asymmetric lattice only. Indeed,
we simulated the original equations of motion (1) with an asymmetry parameter above the critical
value cr (namely  = 0.76) and found that no wave propagation is possible. Using a value slightly
below the bifurcation point (namely  = 0.74) led to a stable wave propagation.
6
Displacement (cm)
0 1 2 3 4
Fo
rc
e 
(N
)
-2
-1
0 
1 
2 
Symmetric variant
R = 22.5, L = 8
r2
r1
(a)
asymmetry parameter ϵ
0 0.2 0.4 0.6 0.8
sp
ee
d 
in 
m
/s
0.4
0.6
0.8
1
1.2
1.4
(b)
Figure 6: (a) Fitted onsite force for the case of R = 22.5 cm and L = 8 cm (dashed red line)
and experimentally measured values (red markers). The symmetric counterpart of this function is
also shown (blue solid line), which was obtained by modifying the local minimum and largest root
along the lines r1 and r2 respectively. The fraction of the distance moved along these lines is the
asymmetry parameter . (b) Plot of the wave speed as the asymmetry parameter  is varied.
Determination of the asymmetry parameter : Let the experimentally measured local minimum
be (u1, F1), the local maximum be (u2, F2), and the largest root be (u3, 0), where the first entry
of each of these coordinate pairs is the displacement and the second entry is the force (see red
markers of Fig. 1a of the main manuscript and of Fig. 6(a)). To explore how the degree of
asymmetry affects wave propagation, we modified the values of x3 and (x2, y2) gradually until a
symmetric function was obtained (see blue curve of Fig. 6(a)). The symmetric function of interest
corresponds to the spline passing through the points (0, 0), (u1, F1), (3u1,−F1) and (4u1, 0) (the
spline was determined with the same procedure used to obtain the solid line of Fig. 1a of the main
manuscript). Let r1 =
√
(u2 − 3u1)2 + (F2 + F1)2 be the distance between the experimentally
measured local minimum (u2, F2) and the value (3u1,−F1). Let r2 =
√
(u3 − 4u1)2 be the distance
between the experimentally measured point (u3, 0) and the value 4u1. We then introduce an
asymmetry parameter  where  = 0 corresponds to the asymmetry of the experimentally measured
onsite-potential and  = 1 corresponds to a symmetric potential. In particular we modify the values
(u2, F2) and u3 such that the distances r1 and r2 change like r1 → r1 − r1 and r2 → r2 − r2 for
various values of  ∈ [0, 1].
References
[1] Izabela K. Kuder, Andres F. Arrieta, and Paolo Ermanni. Design space of embeddable vari-
able stiffness bi-stable elements for morphing applications. Composite Structures, 122:445 –
455, 2015. ISSN 0263-8223. URL http://www.sciencedirect.com/science/article/pii/
S0263822314006424.
[2] H. Xu, P.G. Kevrekidis, and A. Stefanov. Traveling waves and their tails in locally resonant
granular systems. Journal of Physics A: Mathematical and Theoretical, 48(19):195204, 2015.
URL http://stacks.iop.org/1751-8121/48/i=19/a=195204.
7
[3] Thomas Erneux and Gre´goire Nicolis. Propagating waves in discrete bistable reaction-diffusion
systems. Physica D: Nonlinear Phenomena, 67(1-3):237–244, 1993.
8


Unidirectional Transition Waves in Bistable Lattices

https://authors.library.caltech.edu/67867/2/SupplementalMaterial.pdf

Unidirectional Transition Waves in Bistable Lattices

Abstract

Similar works

Full text

Available Versions

Caltech Authors - Main