Effects of Friction Anisotropy on Upward Burrowing Behavior of Soft Robots in Granular Materials

Skins with asymmetric kirigami scales and soft spikes are integrated to the surface of a base self‐burrowing robot, which consists of a soft one‐segment extending actuator. Friction anisotropy is observed at the interfaces between the burrowing robots and different granular materials. Its effects on the pulling resistance and burrowing characteristics are studied. The results demonstrate that the development of friction and friction anisotropy is affected by the characteristics of the granular material, the asymmetric skins, and the relative size of the asymmetric features to the granular particles. Robots with scales or spikes aligned along the upward direction burrow faster than those aligned against the upward direction, especially in relatively coarser granular materials. Particle image velocimetry analysis on the particle displacement fields around the actuator reveals the complexity of dry granular material interactions with soft robots, implying that aligned scales or spikes can impact the distribution of friction preferentially, opening up many possibilities for thoughtful material and geometry‐based manipulation of friction in the design and optimization of future soft burrowing robots for more versatile locomotion capabilities.

Inspired by the burrowing behavior of razor clams, we recently designed and evaluated the upward burrowing behaviors of an extending soft robot in dry sand and developed a soil-mechanics model, which captures the characteristics of the burrowing behavior in shallow (within 7 times the diameter of the actuator) dry sand. [34,38] The robot features a single segment of fiberreinforced silicone tube actuator. The reinforcing fibers limit the actuator to axial extension/contraction motion under pneumatic inflation/deflation. For an actuator vertically buried in sands, cyclic inflation and deflation naturally drive it out of the sand. The natural upward burrowing strategy used by razor clams is intelligent in that it is much simpler than common downward burrowing strategies, which require additional shape-changing features to break the symmetry (as for the dualanchor strategy, peristaltic strategy, or undulation strategy). In this sense, the fact that a one-degree-of-freedom reciprocating motion resulted in net translation itself illustrates passive intelligence which is a consequence of robot-environment interaction. Reflected from the model, the self-burrowing-out behavior is due to the asymmetric soil resistance at each end of the actuator when it undergoes an inflation-deflation cycle.
The model also implies that under the same inflation/ deflation pressure, the burrowing speed of the actuator increases with friction anisotropy, which is defined as the ratio between the upward friction coefficient (resisting downward movement) and the downward friction coefficient (resisting upward movement). This behavior is expected and can be understood intuitively: when the downward friction coefficient is lower than the upward friction efficient, the actuator will experience lower resistance and higher anchorage when it moves upward.
In this article, anisotropic friction features were integrated to the base actuator and their effects on the burrowing behavior were studied. The results highlight the complexity of the interaction between dry GMs with soft robots. We demonstrate that, aside from GM properties, the relative size of the asymmetric surface features (pop-ups and spikes) to the granular particles affect the development of anisotropic friction in granular medium and thus the burrowing behaviors.
The design, fabrication, and characterization of the actuators are first introduced in Section 2.1. The upward burrowing characteristics of three actuators in different GMs and boundary conditions are described in detail in Section 2.2. The results, observed behaviors, and insights from soil mechanics are further discussed in Section 3, followed by the conclusions in Section 4.

Base Actuator
The base actuator is a minimalistic design similar to our previous work, [38] which serves as a reference model for burrowing under a symmetric surface friction feature. The base actuator consists of only one segment of cylindrical silicone tube with double-helix wrapping of Kevlar fibers. The fabrication procedure of the base actuator followed those used elsewhere and a brief description is included in Section 5. [44,45] Multiple actuators were fabricated with one reserved for an experimental control sample; this will be referred to as A1. The remaining actuators were subsequently covered with either kirigami skins or a layer of spiky silicone skins.

Asymmetric Skins
Rafsanjani et al. [39] recently reported that snake-inspired kirigami skins facilitate the crawling of a soft robot on level surfaces. A kirigami skin consists of a plastic sheet with properly designed cuts; when loaded in the designed directions, mechanical instabilities transform the originally flat sheet to an asymmetric 3D-textured surface, which can result in friction anisotropy. Four kinds of cut shapes (e.g., linear, circular, triangular, trapezoidal) were studied in the aforementioned study. It was found that the trapezoidal cuts could propel the soft actuator much more efficiently than the other shapes. In this study, the trapezoidal cut-shape skin design was used and wrapped around the surface of a base actuator. The authors followed the documented design procedures [39] and made adjustments to suit the geometry of the base actuator. In particular, great care was taken to design the cut patterns of the overlapping regions, shown in Figure 1d, to ensure symmetric bonding. Figure 1a shows an enlarged segment of the actuator covered by the kirigami skin, with the scales popped up under inflation. More details on the design and fabrication of the kirigami skins can be found in Section 5 and the Supporting Information.
A thicker layer of the silicone coating used to secure the fibers to the silicone tube was used to embed a layer of asymmetric spikes on the exterior of the actuator. Figure 1c shows the geometry of a single spike. To cover the entire surface of the base actuator, the spikes were distributed evenly along the axial direction with a tip to tip distance of 9 mm and along the radial direction with an inner spacing of about 1.5 mm. Figure 1b shows an enlarged segment of the spiky actuator; more details on the fabrication process can be found in Section 5.
Hereafter, the actuators with kirigami skins are referred to as A21 and A22, where the first numeric digit denotes the actuator type ("2" for kirigami) and the second digit denotes the alignment direction of the actuator ("1" means the scales are aligned against the burrowing direction; and "2" means that the scales are aligned along the direction of the burrowing direction). Two pairs of spiky actuators with different lengths were fabricated, three digits are used to refer the spiky actuators: A311/ A312 and A321/A322, where "31" and "32" differentiate the two spiky actuator pairs and the last digit denotes the alignment of the actuator ("1" for against, and "2" for along). Each pair of experimental tests was based on one actuator oriented in both forward/backward directions; after testing in one direction, the inlet was swapped to the opposite end and tested in the other direction. This treatment ensured that all the actuator characteristics, except the relative direction of the scales/spikes to the upward burrowing direction, were kept constant, isolating the impact of orientation.

Mechanical Properties of the Actuators
Actuation of the actuators was achieved by cyclic opening/ closing the valve using a time-control strategy, [38] which is detailed in the experiment section. The mechanical properties of the actuators were first determined by actuating with periods of 1.8, 3.0, and 3.6 s, whereas the actuators were vertically clamped at one end. Figure 2 shows the deflated and inflated actuators.
The evolution of pressure over time was collected using airpressure sensors located on the control board. The location of the actuator's distal bottom end was captured using a computer vision algorithm based on optical flow. [46] From the pressure and deformation data (see Supporting Information for details), the equivalent Young's modulus was determined by where p is the inflation pressure, ΔL and L are the deformation and the original length, R and r are the outer and inner diameters of the actuator. It was observed that the actuator is hyperelastic and the deformation does not linearly increase with applied pressure. The detailed pressure and deformation curves for the three actuators under three actuated periods are included in the Supporting Information. For simplicity, the equivalent Young's modulus was calculated using the maximum deformation and the corresponding maximum pressure. The results are listed in Table 1. The addition of the kirigami skins greatly increased the modulus of the actuator, due to the high stiffness of the kirigami skin itself. The addition of the spiky layer did not significantly affect the overall modulus, as a similar silicone material was used.

GMs-Actuator Interface Characterization
Three types of GMs were used, namely Ottawa F65 sands (S1 and S1L), Ottawa 20-30 sands (S2), and dry glass beads (GBs). The mean particle size and limit void ratios are summarized in Table 2. S1 was prepared with two relative densities  www.advancedsciencenews.com www.advintellsyst.com (D r ¼ (e max À e)/(e max À e min )). According to geotechnical engineering conventions, [47] the sample with a relative density (D r ) of 40% is considered loose (S1L) and the sample with a relative density of 90% is considered very dense. S2 samples were prepared with relative density of 90% and GB samples were prepared at relative density of 51% (medium dense). The basic characteristics of the GMs are summarized in Table 2. The pullout tests were named based on the actuator type followed by the GM type. The GMs-actuator interface friction was evaluated by pulling out a partially buried dummy actuator, where the motion was induced by external tensile force rather than inner pressure. Using a tensile tester, the total mobilized frictional force was measured as a function of time (F f (t)) at the interface. F f (t) depends on the effective lateral confining pressure ðσ 0 h Þ (or the contact pressure), contact area (A c (t)) and friction coefficient (tan δ), as shown in Equation (2). It is usually convenient to estimate σ 0 h by relating it to the vertical earth pressure σ 0 h via a lateral earth pressure coefficient K, which often varies with initial stress state and the density of the GM. It is often difficult to measure K directly. An alternative approach to characterize friction is to group the effects of K and δ into a single parameter: β ¼ tan δ (the pile design method using β is referred to as "β method" in the field of geotechnical engineering). [48] As the vertical earth pressure σ' v can be easily determined with known effective unit weight γ 0 and depth h ðσ 0 v ¼ γ 0 hÞ, β can be calculated using Equation (2), with the measured total frictional force. A summary of the calculated β is provided in Figure 3 and Table 3.
The pull-out curves for all the GMs-actuator combinations exhibit a similar trend: β initially increases fast to reach a peak value (β peak ) and then gradually drops to a residual value (β residual ) as the travel distance increases. The change in lateral earth pressure (and thus K ) depends on the relative movement of the particles to the contact surface. When the inactive, dummy actuator is pulled out, the GM is subjected to shear loading. If the GM is closely packed (dense), dilation would occur under shearing and the particles close to the dummy actuator will push the surrounding particles located further from the dummy actuator. In this process, the normal contact pressure acting on the dummy actuator surface will increase. Another way to look at it is to assume the surrounding soil behaves like a spring with one end fixed at the far field and the other connected to a small soil element in contact with the dummy actuator. When this soil element dilates,    it pushes the spring away from the dummy actuator and thus increases its loading on the robot's surface, assuming the stiffness of the spring is a constant. After a certain amount of movement, the soil close to the dummy actuator yields and the friction reaches its limiting value. In contrast, if the soil is loose and hence contractive, particles close to the actuator would be compressed during shearing, and the soil element would move toward the dummy actuator surface, the spring would be unloaded, resulting in decreased contact pressure. The significant difference between β peak and β residual observed in the denser materials reflects the dilative nature of the materials. Furthermore, for a dummy actuator with finite length, when it is pulled upward, voids created around and below (i.e., bottom end) it will allow soil to flow down due to gravity, filling the voids created. This causes further pressure relief (decreasing K and σ 0 h ) on the dummy actuator. This process complicates the friction phenomena and differentiates the GM-structure interaction from classical stick-slip processes observed between rigid surfaces. The vertical displacement required to reach β peak is small and typically smaller than 1 mm (Figure 3a,b). The results show that the interface friction is affected by the density of the GM, the particle size and shape, as well as surface texture and the direction of the asymmetric features. The increase in relative density of S1 from 40% to 90% resulted in threefold increase in β peak due to the less dilative nature of the looser material; this is also consistent with the fact that the lateral earth pressure coefficient K increases with relative density. [49] For GMs with the same relative density, experiments in Ottawa sand showed greater β peak values than in glass beads, which can be attributed to the smooth, spherical shape of the glass beads. Larger particle sizes led to higher frictional anisotropy (β against /β along ). The spiky actuators exhibit highest anisotropy of β residual in glass beads (%1 mm), followed by Ottawa 20-30 (%0.72 mm). Pulling the spiky actuators out from Ottawa F65 (%0.20 mm) did not show obvious anisotropy.

Burrowing Tests
A series of burrowing tests were conducted in a cylindrical container to evaluate the burrowing performance of actuators with different surface features. The experiment setup used for the 3D burrowing tests is shown in Figure 4a. A metal tube was connected to a vertically aligned soft actuator, which was pre-buried inside the soils. A pneumatic control board regulated the opening and closing of the valve, leading to the inflation/ deflation of the actuator, which burrowed up out of the soil gradually. Two cameras were used: one to record the soil surface deformation during the test, and another to monitor the burrowing progress by tracking markers on the metal tube. Representative frames from the sequence are shown in Figure 4b. All the burrowing tests followed a standardized procedure, which was repeated three times. For each burrowing scenario, the soils were prepared with consistent relative densities and a same actuation period of 3.6 s. All the results for the burrowing-out of each actuator showed great consistency and repeatability. Details regarding the standard testing procedure can be found in Section 5 and the Supporting Information.
A similar 2D setup was used to capture the soil displacement field associated with the burrowing actuator in a plane-strain condition, as shown in Figure 4c. The soil displacement corresponding to the actuator inflation/deflation during burrowing is videoed. The particle image velocimetry (PIV) technique was implemented to recover the soil displacement field from the recorded image flow. Details regarding the 2D test procedure and the PIV analysis are provided in Section 5. An example video can be found in the Supporting Information. Figure 5a shows a complete burrowing curve, which exhibits a general "S" shape and includes multiple burrowing cycles. Each cycle of burrowing can be considered a single gait cycle, in which the top end of the actuator first moves up (advances) when inflated and moves slightly down (slips) when deflated. The corresponding travel distances are referred to as advancements and slips, respectively; and the difference between each pair of advancement and slip is considered as the stride length. Figure 5b shows the characteristics of the burrowing process of  www.advancedsciencenews.com www.advintellsyst.com the case A1S1. The rates of advancement and slip per cycle both increase over a number of cycles each trial, which can be attributed to the decreasing resistance required for the actuator to overcome the decreasing embedment depths as it moves upward. The advancement increases faster during initial cycles and slip increases faster near the end of a trail. This causes the stride length per cycle to increase first and then decrease.

Base Actuator (A1S1)
Comparing the stride length curve to the burrowing curve and analyzing the top-view videos, it was found that the turning point of the stride length curve occurred at the 12th cycle, 5 cycles before the top of the actuator emerged from the sand.

Kirigami Skinned Actuator (A21S1/A22S1)
The general shape of the burrowing curves and characteristic curves for the actuators covered with kirigami skins (A21 and A22) are similar to that of A1S1. A21 and A22 burrowed much more slowly in S1 than A1 ( Figure 6), which can be attributed to the much higher Young's modulus and the higher frictional resistance induced by the pop-up features of the skin. Friction anisotropy resulted in higher stride lengths and thus higher burrowing speed for A22, which has scales aligned with the upward direction. The anisotropy effect did not appear until the 12th cycle, before which the advancements and slips for actuators A21 and A22 are almost identical to each other. The authors believe that at the beginning stage, the axial deformation of the actuators was small (lower than 4 mm), preventing the kirigami skin from expanding more than about 0.4 mm, resulting in a low anisotropy. Afterward, the scales expanded more, and the corresponding higher friction resistance led to lower advancements and slips, but the stride lengths kept increasing due to decreasing soil pressure. The anisotropy caused higher stride lengths for A22, which took about 90 s (or 25 cycles less) to burrow out of the sand.
Additional reasoning for this behavior is included later in Section 2.4. The burrowing characteristics of spiky actuator (a) (b) Figure 6. a) Burrowing curves for cases A21S1 and A22S1. b) Burrowing characteristics for cases A21S1 and A22S1. The data points reflect the average value of the characteristics, whereas the shaded area covers 1 AE standard deviation.
www.advancedsciencenews.com www.advintellsyst.com pairs (A311S1L/A312S1L, A321S1/A322S1) did not show noticeable effects from the orientation of the spikes, which was consistent with the low friction anisotropy of spiky surfaces in S1, as observed in the pull-out tests ( Figure 3c and Table 3). The only notable differences were found in the first 6 cycles for tests in loose Ottawa F65 (S1L), where actuator A22 experienced smaller slips than A21. Note that after the 6th cycle, the top end of A22 reached the soil surface.
Ottawa 20-30 (A321S2 and A322S2): The burrowing characteristics of the spiky actuators in soil S2 (medium dense Ottawa 20-30) were significantly different from those in S1. There were "breaking-through" moments, before which the stride lengths were extremely small (%0.2 mm) and after which the stride lengths significantly increased ( Figure 8). It will be shown later in Section 2.4.1. that the breaking-through behavior is related to soil failure surfaces. The "breaking-through" Figure 7. a) Burrowing curves for cases A311S1L and A312S1L; b) Burrowing characteristics for cases A311S1L and A312S1L. c) Burrowing curves for cases A321S1and A322S1; d) Burrowing characteristics for cases A321S1and A322S1. The data points reflect the average value of the characteristics, whereas the shaded area covers 1 AE standard deviation.
(a) (b) Figure 8. a) Burrowing curves for cases A321S2 and A322S2; b) Burrowing characteristics for cases A321S2 and A322S2. The data points reflect the average value of the characteristics, whereas the shaded area covers 1 AE standard deviation.
www.advancedsciencenews.com www.advintellsyst.com moments for the actuator A322 (at 198 s and 95 mm deep) occurred earlier than A321 (at 227 s and 88.5 mm deep); after breaking through, A322 burrowed faster with larger stride lengths than A321. The earlier breaking-through moment and the larger stride lengths afterward indicated the effects of friction anisotropy: the friction force experienced by A322 was lower than that of A321 in medium dense Ottawa 20-30 sands. This trend was also consistent with the higher friction anisotropy of spiky surfaces in S2, as observed in the pull-out tests ( Figure 3 and Table 3). Based on the comparison between the behaviors of spiky actuators in S1 and S2, it was hypothesized that the particle size affects the development of friction anisotropy and thus the burrowing behavior. To validate the hypothesis, burrowing tests were conducted in glass beads, which have even larger particle size. To further understand the effects of friction anisotropy on the burrowing behaviors of the spiky actuators, experiments were conducted at plane-strain conditions to visualize the particle movements around the actuators.
Plane-strain tests in Ottawa 20-30 (A321S2-2D and A322S2-2D): The plane-strain condition burrowing curves of the spiky actuators in Ottawa 20-30 sands (Figure 9) resemble that in the 3D cylindrical container (Figure 8). In the 2D chamber, the actuators were embedded deeper (about 320 mm) and this led to longer total burrowing time. The breaking-through moments for both actuators occurred at a similar depth (about 123 mm), which was deeper than the 3D cases. The stride lengths for A322 were greater than those for A321 at the same embedment depths, which was due to greater advancements and smaller slips for A322. This led to three times faster burrowing speed for A322. These trends confirmed the effects of anisotropic friction.
Similar trends were also observed for the tests in dry glass beads ( Figure 10). The effects of anisotropic friction were more significant, indicated by greater differences in stride lengths between A322 and A321. This is also consistent with the highest friction anisotropy of spiky surfaces in glass beads, as observed in the pull-out tests (Figure 3 and Table 3); and it further proved that larger particle size enhanced the development of friction anisotropy. The breaking-through moments occurred at a deeper depth (about 140 mm), probably due to the overall lower friction in glass beads and the lower specific gravity of glass compared with that of quartz. Before breaking through, A322 burrowed faster in glass beads than in Ottawa 20-30, but it burrowed more slowly afterward. This is due to the fact that lower friction not only increases the advancements but also results in increase in slips (Figure 9 and 10). www.advancedsciencenews.com www.advintellsyst.com

Particle Displacement Fields
In general, the development of anisotropic friction during the burrowing processes agreed with that observed in the pull-out tests. As discussed in Section 2.1.3., the mobilized friction between a moving object and the surrounding GM depends on the dilatancy of the material or the movements of the particles close to the object. Particle image velocimetry (PIV) method was used to retrieve the displacement fields of the particles interacting with the upward burrowing robot. The PIV technique has been widely used for the measurement of soil deformation in the field of geotechnical engineering. [21,[50][51][52][53][54][55] It recovers the soil displacement field between image pairs by dividing the initial image into a mesh of square subsets and tracking the displacement of image subsets by searching the maximum correlation between the initial and target subsets. In this study, the PIV software GeoPIV-RG was implemented for the soil displacement field analysis. [56] GeoPIV-RG adopted a reliability guided computation method to enhance the accuracy and efficiency of small and large soil deformation measurement. Details regarding the GeoPIV-RG software can be found elsewhere and the parameters used in this study can be found in Section 5. [56] Four individual burrowing cycles were selected to analyze and compare the evolution of the displacement fields. These included two cycles from test A321S2-2D. One cycle corresponded to the one started at the 137th second, representing burrowing at a deep location (embedment depth of about 260 mm). The other cycle corresponded to the one started at the 1070th second, representing the cycle at the breaking-through moment (embedment depth of about 125 mm). Two comparable cycles were selected from test A322S2-2D, corresponding to cycles started at the 74th and 363rd second. Figure 11 shows the accumulated vertical displacement contours during the inflation phase of the four selected cycles. The top and bottom ends of the actuator roughly locate at the points with minimum and maximum vertical displacement, respectively. For all the four cycles, the vertical displacement fields feature two distinctive regions: an upward moving region around the top segment of the actuator, and a downward moving region around the lower segment of the actuator. Between these two regions, there are thin segments with negligible movements. These intermediate segments represent the anchor points during inflation. Above the anchor points, the actuators move upward, dragging or pushing the surrounding particles upward. Below the anchor points, the opposite occurs. In general, actuator A322 (with spikes aligned along the upward direction) induced higher vertical displacement magnitudes, indicating greater extension than A321. The displacement fields for the deeper burrowing cycles are more localized, with bulb-shaped influence zones at the two ends of the actuator. The influence zone above the actuator extended to about two times of the diameter of the actuator, which agrees with the observations on pulling out deeply embedded plate anchors from loose sands. [52] The influence zone at the bottom end of the actuators extended to one diameter of the actuator and resembled that of a penetrating cone or pile. [57] The influence zone from the breaking-through cycles extended to the soil surface and forms a truncated cone, analogous to that of a shallow plate anchor being pulled out. [52] The horizontal displacement fields ( Figure 12) for the inflation phases have four distinctive regions. From bottom up, the first region is located below the bottom end of the actuator, where the soil particles are pushed sideways. The second region is located between the bottom end and the anchor point, where the soil particles move toward the actuator, indicative of decreasing contact pressure. The third region is located between the anchor point and the top end of the actuator, where the soil particles moves away from the actuator, indicative of increasing contact pressure. Above the top end is the fourth region, where the soil is pushed sideways significantly. In general, the inflation of A322 resulted in stronger horizontal movements of the surrounding particles, indicating more significant change in contact pressures. www.advancedsciencenews.com www.advintellsyst.com

Deflation Processes
Caution is needed when interpreting the accumulated vertical displacement fields for the deflation processes ( Figure 13). Unlike the inflation processes, the two ends of the actuator and the anchor points could not be identified easily from the displacement fields. The downward moving region for all cycles extended to the soil surface. The upward moving region may not be captured well, as in the breaking-through cycle of A321. When it is indeed captured, it may indicate the movement of the bottom end of the actuator or the "void" space beneath the bottom end. The horizontal displacement fields ( Figure 14) indicate that during the deflation phase the soil particles around the actuators move toward and alongside the entire actuator body, implying decreasing contact pressure and thus lower friction resistance. Again, the A322 induced stronger motion of particles than A321.
The displacements of the top end of the actuators (d t ) for the selected four cycles can be obtained from the burrowing curves (Figure 9a). Similarly, the displacements of the bottom end of the actuators (d b ) were obtained using the optical flow algorithm. The results are summarized in Table 4.
From the PIV analysis, it is clear that the locomotion of the robot greatly affects the surrounding soil. The loading and unloading of the soil cause changes in soil stress and packing states; consequently, the interface frictional behavior is affected. [54][55][56][57][58][59][60][61] All observations so far confirmed that the friction anisotropy can develop in GMs and such anisotropy affects the burrowing characteristics of the robots. In Section 3, a more indepth discussion on the soil-actuator friction is provided, which shed light on the development of future self-burrowing robots.  www.advancedsciencenews.com www.advintellsyst.com 3. Discussions

Soil-Robot Interface Friction
With the contact pressure changing on the robot surface during the burrowing process, it is difficult to directly quantify the soil-robot interfacial friction. Here we attempt to back calculate the equivalent β value (and thus the mobilized friction) during the inflation phases based on soil mechanics and the observed kinematics of the actuators. The general idea of the proposed approach is to simplify the inflation phase as a quasi-static process, where all the external forces on the actuator are balanced (Equation (3)) The external forces experienced by the actuator include the forces from the sand at the top and bottom of the actuator, F t and F b , respectively, the gravitational force G, and the downward and upward frictional force from the sand, F df and F uf .
The tensile force along the actuator ( Figure 15) is then where x is the distance from the soil surface to the point of interest, A in is the inner cross-sectional area of the actuator, P in is the inflation pressure, L is the original length of the actuator, F f (x) is the total frictional force from the top end of the actuator to the point of interest, and h is the embedment depth measured from the soil surface to the top end of the actuator. Assume that the distance from the top end of the actuator to the anchor point is x a . At locations above the anchor point At locations below the anchor point  Positive values indicate upward movement; b) measured from the bottom end of the actuator.
www.advancedsciencenews.com www.advintellsyst.com The total movement of the top end (d t ) and bottom end (d b ) of the actuator can then be calculated as TðxÞ EA a dx (7) where E is the equivalent modulus of the actuator and A a is the actuator wall area. The values of h, d t , d b, and x a are given in Table 4; L, E, A a , A in are given or can be calculated based on Table 1; β 0 can be calculated based on data in Table 2; P in is 207 kPa, and G is 0.5 N. To back calculate β 1 and β 2, an estimation of F t should be made. Based on the PIV analysis, it was found that the displacement field around the top end of the actuator during the inflation phase resembles that of a plate anchor being pulled upward. Therefore, F t is analogous to the pull-out capacity of an embedded anchor. In geotechnical engineering, the pull-out capacity of an anchor embedded in soil are usually calculated using a breakout factor N γ , which is related to the soil strength parameters, anchor geometry, and its embedment depth. [49] Thus, it is reasonable to use the breakout factor N γ to estimate F t , or, where A is the total cross-sectional area of the actuator. In this study, the expression used is derived for pull-out capacity of plate anchors buried in sand under plane-strain conditions. [49] To calculate N γ , the dilation angle ψ of Ottawa sands is estimated using the dilatancy relationship, [62] which is a modified version of the original relationship to account for low confining pressure levels. [63] More details regarding the calculation of N γ can be found in the Supporting Information.
The back calculated values of β for the inflation processes of the selected four cycles are listed in Table 5, along with the estimated β values based on the pull-out tests. The estimated β values are calculated by averaging the direction-dependent β values, ranging from zero to the deformation of the upper or lower segments of the actuator. For example, β 1 for the inflation phase of the cycle A321-191 s is calculated by averaging the pull-out curve (A321S2 in Figure 3c) from 0 to 3.64 mm, which is the movement of the top end of the actuator d t ¼ 3.64 mm (see Table 4); to calculate β 2 , the pull-out curve for A322S2 is used, as the lower segment of the actuator moves downward, and the spikes align along the moving direction.
From Table 5, it is clear that higher β values were developed along the actuator A321 than along A322, indicative of friction anisotropy. In addition, higher β values were developed at shallower depths, which is consistent with the fact that the lateral earth pressure coefficient K increases with decreasing confining pressure. [49,64] At greater depths, the back calculated β 1 values are greater than those of the estimated values based on the pull-out data, whereas they are lower than the estimated values at shallower depths. It is likely that the robot kinematics alter at different depths, due to the varying pressure build-up around the actuators. In contrast, the back calculated β 2 values are much lower than those of the average values estimated based on the pull-out data. This is due to the movements of the soil particles toward the lower segments of the actuator, which greatly reduces the contact pressure.
The aforementioned analysis confirms the complex nature of friction and friction anisotropy around a soft burrowing robot. The mobilized friction directly depends on the kinematics of the robot and on the resulting particle movements.

(Relative) Particle Size Effects on Friction Anisotropy
For spiky actuators, friction anisotropy was not observed in the pull-out tests or the burrowing tests in S1 (Ottawa F65), but was observed in S2 (Ottawa 20-30) and glass beads (GBs). The mean particle sizes for the three materials are 0.2, 0.72, and 1.0 mm. The results suggest that friction anisotropy increases with particle size. A more proper hypothesis may be developed on the relative size of asymmetric surface features to particle, instead of the absolute sizes. In fact, the actuator covered with kirigami skins also showed clear friction anisotropy in the finer Ottawa F65. Note that the geometry of the kirigami scales changes with the strain level. The height of the scales increases from a minimum to maximum value during the whole course of inflation; it then decreases during deflation. Ideally, the kirigami  Figure 15. Schematic for the model used to back calculate the β values. a) free body diagram of the actuator; b) free body diagram on part of the actuator.
www.advancedsciencenews.com www.advintellsyst.com scales will fully close when deflated. However, it was observed that sand particles could be trapped inside, leading to jamming. In addition, since the strain was not developed uniformly along the entire actuator during inflation, the height of the scales would vary from location to location. Thus, an accurate estimation of the scales' geometrical characteristics could not be made. Nevertheless, it was observed that the maximum height of the kirigami scales was 1.37 mm under inflation pressure of 207 kPa with no confinement; this corresponds to 6.85 times of the mean size of Ottawa F65. The real height of the scales in sand is expected to be much smaller than this limit value, due to the small deformations and confinement from the soil. In contrast, the spikes had a constant height of 2.5 mm and were molded directly onto the base actuator. During inflation and deflation, the change in height of the spikes was negligible, but the spacing of spikes along the longitudinal direction may change due to the extension and contraction of the actuator. The height of the spikes were 12.5 times of the mean size of Ottawa F65, 3.5 times of the mean size of Ottawa 20-30, and 2.5 times of the mean size of glass beads. Considering all the cases, there may exist a threshold spike-to-particle ratio, between 6.85 to 12.5, to induce friction anisotropy. But this assessment should be viewed with caution, due to the limited number and variations of the scale sizes tested in this study. Size effects on the development of friction anisotropy at soil-structure interfaces were also hinted in another recent study. [22,65] Interface shear-box tests were conducted to characterize the shear resistances (friction) at the interface of soil and snakeskin-inspired surfaces. The effects of the "scales" on friction and friction anisotropy were systematically studied by varying the length and height of the "scales" on the surface. The scale height used in their study varied from 0.1 to 0.72 mm, much smaller than that used in this study; the scale lengths varied between 6 and 32 mm. Their testing results suggest that the residual friction anisotropy for scales with height of 0.3 and 0.7 mm developed in Ottawa 20-30 was higher than that developed in Ottawa F65. It was shown that although the mobilized friction in both sands decreased with increasing scale lengths from 6 to 32 mm (while keeping the height at a constant of 0.3 mm), the friction anisotropy in Ottawa 20-30 was not significantly affected. The friction anisotropy in Ottawa F65 decreased with the scale length and eventually diminished when the scale length was 32 mm. In this study, the length of the kirigami pop-up scale was about 2.0 mm and that of the spikes were 7.5 mm, well within the range where friction anisotropy was observed in the previous study. The results from the two studies cannot be compared directly, as the asymmetric surface features in the two studies have different geometries, patterns, and compliances. However, both suggest the effect of relative scale-to-particle size may play a role in the development of friction anisotropy. It warrants further in-depth study on this matter in the future.

Implications and Outlook
The results of the pull-out tests and burrowing tests reported herein demonstrated the effect of asymmetric surface features on the mobilized friction in GMs. The robots used in this study burrow upward as a result of the relatively lower overall resistance in the downward direction. However, the findings from this study can shed light on the design of future soft burrowing robots, which can burrow in other directions (e.g., downward). For example, to burrow downward, the overall resistance in the downward direction should be higher than in the upward direction. A potential approach is to induce higher side friction resistance in the downward direction with spikes/scales aligned against the upward direction (or along the downward direction). This requires a more systematic study on friction and friction anisotropy in GMs.
In most cases, introducing roughness to an originally smooth surface will result in higher friction resistance, no matter if the roughness features are symmetric or asymmetric (compare A3 cases to the A1 case in Table 2). Note that the roughened surface has the same chemical composition as the original surface. The kirigami skin used in this study, however, resulted in smaller peak β values compared with A1. This is due to the fact that the undeformed plastic sheet used in this study is very smooth, and induces lower friction than the silicon surface. Also, the scale-like pop-out feature in current design aligns with the longitudinal direction of the actuator and tends to compress the sand behind during inflation; the pop-up may be constrained because of the stiffness difference between the sands and the skin material. Therefore, the absolute magnitude of the friction values can be controlled by modifying the skin materials and the orientation and alignment of the pop-out.
In addition, the stiffness of the asymmetric elements as well as the joints connecting them to the main body affects the development of friction and friction anisotropy. [20] Friction anisotropy inversion, or the development of negative friction anisotropy, may occur when the dominating mechanism switch from mechanical interlocking to adhesion-mediated friction, or vice versa. One example is that high-aspect-ratio soft scales exhibits negative friction anisotropy when sheared along/against a soft substrate. GM may be seen as a "soft" deformable surface with different levels of roughness depending on particle size distribution and density. This implies inversion of friction anisotropy is possible in GMs. In fact, when the kirigami dummy actuator was pulled out against the scales demonstrated a slightly negative anisotropy, where a lower residual β value was measured (Figure 3b).
Collectively, interface friction could be controlled to modify the absolute magnitude of the friction and degree of anisotropy, as well as to induce reversion of anisotropy as needed. Recent studies on control of frictional asperities using magnetic or electric means opened new opportunities to manipulate friction properties actively. [66] By leveraging on friction manipulation and embodying material intelligence, endless potentials exist to advance the future designs of soft burrowing robots for more versatile locomotors in GMs.

Conclusion
Skins with asymmetrical surface features (kirigami pop-up scales and spikes) were used to cover a soft burrowing robot, in the hope of introducing frictional anisotropy at the interface between the robot and GMs. The results showed that the development of friction anisotropy at the interfaces were affected by the characteristics of both the surface geometry and the GM. Notably, anisotropy developed at the interfaces between the kirigami skinned actuators and the finer sand, at the interfaces between the spiky actuators and the coarser GMs, but not at the interface between the spiky actuators and the finer sand. The developed friction anisotropy affected the upward burrowing characteristics: robots with scales or spikes aligned along the upward direction burrowed faster than those aligned against the upward direction. When deeply buried, the actuators initially burrow slowly until it reaches the breaking-through depth (about 7 times of the diameter), after which the burrowing speed significantly increases. Using PIV analysis, the sand displacement fields induced by the burrowing robot revealed complex but consistent particle flow patterns around the actuator. During inflation, the materials located close to the two ends of the actuator and along the side walls above the anchor point were pushed away, indicating pressure build-up and thus increased downward friction resistance. The materials located along the robot surface below the anchor point moves toward the robot, indicating pressure relief and thus decreased upward friction resistance. Upon deflation, particles flow mostly downward around the actuator and would fill in the void below created by the upward contraction of the lower segment of the actuator. The PIV analysis also revealed that before breaking through, the particle movement was more localized. However, at the breaking-through movement, the displacement zone extended to the soil surface. These trends are reflected by the back calculated values of parameter β, which combines the effect of lateral earth pressure coefficient and the friction coefficient. Incorporating friction anisotropy in the design of future soft robots is potentially a promising approach to generate anchorage needed for burrowing in other directions.

Experimental Section
Fabrication of the Based Actuators: In brief, the two-part liquid silicone Dragon Skin 30 (Smooth-On) were selected for the molding of the main body of the actuator. A symmetrical, double-helix wrapping of Kevlar fibers at angles of AE 75 to the longitudinal axis of the actuator was then applied to the silicon tube, ensuring a pure extension/contraction at the pressurization/depressurization. Afterward a thin layer of liquid silicone (1.45 mm thick, Dragon Skin 10, Smooth-On) was applied to the wrapped silicone tube to secure the threads. Once cured, both open ends of the silicon tube were sealed using liquid silicone, and a quick-release tube-fitting adaptor was inserted into the actuator serving as an air inlet/outlet. Silicone glue was applied to secure the adaptor and to avoid air leakage.
Fabrication of the Kirigami Actuators: Care was taken when designing the cut patterns to ensure that the kirigami skin can be attached to the base actuator without creating asymmetrical deformation in the axial direction. It was decided to use two separate skins with the same width to cover the actuator. As shown in Figure S1, Supporting Information, each half of the skin was divided into three zones with different cut patterns. The cut pattern shown in Figure S1a, Supporting Information, was used for the central region, which will be in contact with the base actuator. The cut pattern shown in Figure S1b, Supporting Information, was used for the side regions to facilitate assembling. The cut pattern used in the side regions have a greater spacing (2.4 mm) along the direction of 60 to the horizontal, than that in the central region (1.0 mm). This design ensured that the overall stiffness of the kirigami skin was not increased by the side regions when overlapped. To attach the kirigami skins to the base actuator, the procedure shown in the Supporting Information was followed.
Fabrication of the Spiky Actuators: After wrapping the inner silicone tube with fibers, a thin layer of liquid silicone (%1 mm thick, Dragon Skin 30, Smooth-On) was first applied. The resulting tube was then immediately slid into a 3D-printed mold with negative spiky patterns on the inner surface. After the thin layer of liquid silicone became squishy (%2 h after silicon application), more vacuumed liquid silicone (Dragon Skin 30, Smooth-On) was poured into the mold. This method was proved by the authors to create a firm bonding between the inner silicone tube and the external skin layer. Other fabrication steps were similar to those used in fabricating the base actuator.
Robotic Control: An open-source pneumatic control board was built to drive the actuator in this study. [67] The control board featured a microcontroller (Arduino Mega), miniature diaphragm air pump, 4-channel metal-oxide-semiconductor field-effect transistor (MOSFET), multiple pneumatic solenoid valves, and pressure sensors. For simplicity, a time-control strategy was adopted for the actuator alternation behavior. Specifically, we controlled the opening and closing time of the valves to achieve inflation (extending) and deflation (contracting) of the actuator. Note that the pump was running at its rated voltage, so that the pressurized air source did not change with time. During the actuating period T, the valve was opened for the first one-third of the time and closed for the remaining two-third of T. For the calibration tests in the Supporting Information, three periods (1.8, 3.0, and 3.6 s) were used to inflate the actuator with different pressure levels (150, 187, and 207 kPa, correspondingly). For all the burrowing tests, a period of 3.6 s was used, resulting in a maximum pressure of 207 kPa.
Pull-Out Tests: A tensile strength testing device (MTS E42.503 Crosshead) was used to conduct the pull-out tests. Three to five tests were conducted for each scenario and the average values and standard deviations were reported in Figure 3 and Table 2. The detailed standard procedure is further described in the Supporting Information.
3D Burrowing Tests: The burrowing tests were conducted in a cylindrical bucket which measures 430 mm in height and 406.4 mm in diameter. The tests also followed a standard procedure to ensure repeatability. A step-bystep description can be found in the Supporting Information. In brief, the soil samples were first prepared by a dry pluviation (sand raining) method to achieve a relatively loose state. A vibrator (VP51D1) was used to densify the soil to the desired relative density by monitoring the surface level of the soil. Due to the lack of surcharge loading on the top surface during vibration, the prepared samples were not homogeneous in the vertical direction. Specifically, the void ratio decreased with depth slightly (within 6%, see Supporting Information for a more detailed quantitative study on the variation of void ratios). Nevertheless, the samples were prepared in a consistent manner and the results from different trials are comparable. During the burrowing-out tests, videos were taken from the side and top of the container using two DSLR cameras (Nikon, D3400) at 1920 Â 1080 pixel and 24 fps. The detailed vertical moving trace of the actuator was then obtained from the side-view video using an open-source computer vision library OpenCV. [46] An optical flow algorithm based on the Lucas-Kanade method was used to track the motion of a marker located on the metal tube. [68] The scaling factor, which relates the distance measured by pixel using the optical flow algorithm to the physical length, was determined through eight evenly distributed markers on the metal tube. For each burrowing scenario, three burrowing tests were conducted to ensure the repeatability of the results.
2D Burrowing Tests: The 2D setup was similar to that used for the 3D burrowing-out test, except that a transparent thin chamber was used to model the plane-strain condition (see Figure 4c). The preparation and testing procedure for the 2D test followed that of the cylindrical container tests. Each frame of the video for 2D tests were processed using PIV technique, from which the entire displacement field was obtained.
PIV Analysis: In this study, the soil deformation analysis focused on the local area around the burrowing robot. All the analyses were performed under a minimum reference-and-target subset correlation coefficient of 0.9 and a minimum full-field correlation of 0.75. Due to the relatively large involved deformation, the adopted mesh sizes varied between 50 and 75 pixels, corresponding to 18-27 mm.