Imaging dynamic three-dimensional traction stresses.

Traction stress between contact objects is ubiquitous and crucial for various physical, biological, and engineering processes such as momentum transfer, tactile perception, and mechanical reliability. Newly developed techniques including electronic skin or traction force microscopy enable traction stress measurement. However, measuring the three-dimensional distribution during a dynamic process remains challenging. Here, we demonstrated a method based on stereo vision to measure three-dimensional traction stress with high spatial and temporal resolution. It showed the ability to image the two-stage adhesion failure of bionic microarrays and display the contribution of elastic resistance and adhesive traction to rolling friction at different contact regions. It also revealed the distributed sucking and sealing effect of the concavity pedal waves that propelled a snail crawling in the horizontal, vertical, and upside-down directions. We expected that the method would advance the understanding of various interfacial phenomena and greatly benefit related applications across physics, biology, and robotics.

INTRODUCTION

Traction or contact stress between objects plays an important role in nature and industry including earthquake caused by crustal movement (), locomotion of organism (), human dexterous operation (), and performance of bearings and gears (). Obtaining the distribution and evolution of interfacial traction stresses is the basis for understanding the natural phenomena and improving engineering performances. For the last hundreds of years, numerous theoretical models have been proposed to quantify the magnitude and distribution of traction stresses with considering pure elasticity (), adhesion (, ), friction (), surface roughness (, ), and viscoelasticity and plasticity (, ). Recently, by means of the development of computer simulations, the theoretical framework for solving contact problems has been notably improved (, ). However, the progress of traction stress measurement technology is limited but highly demanded.

One of the typical examples to detect the temporal and spatial information of traction stresses is the human finger, which therefore can precisely recognize surface texture () and reliably grab objects (). Recently, functional devices that mimic properties of human skin known as electronic skins have been greatly developed (). In general, the pressure distribution can be easily obtained through arrays of piezoresistive () or capacitive sensor () or qualitatively represented by mechanoluminescence materials () and photoelastic materials (). By means of complicated electrode array structures, both normal and shear traction forces are available to be measured with a certain spatial resolution (–). These electronic skins with other potential notable abilities such as flexibility, self-healing, and wireless communication have great application prospects in wearable devices, robotics, and prosthetics (, ). However, densely packed sensor arrays may cause difficulties in signal processing, wiring fabrication, force calibration, and achieving high spatial resolution ().

One alternative approach to measuring the distribution of traction force is based on optical method, which allows us to convert it into an image and take the advantage of the updates of high-resolution cameras and image-processing algorithms (). One of state-of-the-art method that can infer the shear stress distribution was called traction force microscopy and commonly used in the mechanics of biological cells (, ). It is generally achieved by tracking the in-plane movement of fluorescent particles embedded in the compliant substrate, which follows a similar principle to the particle image velocimetry technique for flow visualization (). An improvement to obtain the three-dimensional (3D) traction forces can be achieved by layer-by-layer scanning using confocal microscope (, ), which, however, was only applicable to small deformation (approximately micrometers) and slow variation (approximately seconds to minutes). Thus, the application becomes limited when involving large deformation and dynamic process.

Facing the great demand and challenges for characterizing 3D interfacial traction stresses with both high temporal and spatial resolution, we further extended the traction force microscopy technique based on the binocular stereo vision technology. By means of the specially designed multilayer layout, a binocular camera and advanced global feature matching algorithm called digital image correlation (DIC) were used to capture and reconstruct the 3D deformation of the contact surface. The three-axis traction stresses with a typical spatial resolution of about 10 μm and a temporal resolution of about 10 ms can be then obtained by deconvolution from the deformation field based on mechanical models. We introduced an iteration method to perform the deconvolution with high numerical stability. The resolution and field of view can be readily adjusted by adjusting the optical system as needed. The proposed cost-effective method was demonstrated to be powerful for understanding various interfacial phenomena and had the potential to become a tactile sensor for using in robotics. For example, we used it to image the two-stage adhesion failure process of microarrays and displayed the contribution of both elastic resistance and adhesive traction to rolling friction at different contact regions. Furthermore, it successfully revealed the crawling mechanism of a snail that the distributed sucking and sealing effect with the roll-in and peel-off behaviors of the pedal waves propelled a snail crawling in the horizontal, vertical, and upside-down directions. The proposed method has unique advantages in providing traction stress information of contact interface with high spatial and temporal resolution, which can have great application prospects in various scientific and industrial fields.

RESULTS

Principle of traction force measurement

The basic principle of the method is similar to the traction force microscopy but with several differences to meet the demands of dynamic 3D stress measurement. One of the most important one is to use binocular cameras providing stereo imaging capabilities, which can overcome the shortage of slow scanning speed of a confocal microscope. It also has the potential to become a tactile force sensor for using in robotics after the miniaturization and improvement of real-time performance. Moreover, we introduced an iterative method to calculate the traction stress from measured deformation under the surface, which can avoid directly calculating the inverse of ill-conditioned matrix. It can substantially improve the numerical stability, avoid high-frequency errors amplified by the inversion process, and reduce the dependence of the results on low-frequency pass filter.

The structure layout is designed with several layers (shown in Fig. 1A). A thin reflective layer covering on the top surface of a transparent elastomer block (silicone rubber; Young’s modulus of 0.04 and 0.095 MPa were used in this work) to create a clean internal light environment. The second layer closely under the reflective layer is the dense speckle layer. It can provide the matching features to reconstruct the shape and deformation. The elastomer block was bonded to a glass plate, and commercial binocular cameras were fixed below the plate, which can capture image frames up to 90 frames/s (fps) at a resolution of 640 pixels by 480 pixels. The two cameras were fixed at an angle of about 15° from the vertical direction as a trade-off between the overlap of view and the depth of field.

Fig. 1.

Outline of the dynamic traction stress imaging method.

(A) Schematic illustration of the multilayer structure. It can identify the ratchet-like microstructure from shear stress distribution when pulled in different directions. (B) Principle of obtaining the interfacial stress from the captured figures of the deformed surface.

When the traction force is applied on the free surface of the elastomer, the force information is transformed into surface deformation. The feature layer under the top surface can be used to reconstruct the deformation (shown in Fig. 1B). One of the state-of-the-art methods to achieve this target is the 3D-DIC algorithm (). Theoretically, if the shape, boundary condition, and constitutive equation of an elastic body were predetermined, then the relationship between deformation and applied force can be uniquely determined. In practice, one reasonable simplification to apply Boussinesq’s solution of elastostatic equation is by assuming that the elastic body is homogeneous and semi-infinite (SI) when the deformation and contact scale are both smaller than the thickness (). In this case, the deformation field u (i, j = x, y, z) takes the form of u(z0) = K(z0) * p with summation over repeated indices, where the symbol “*” indicated convolution, z0 is the depth of the speckle feature layer from the surface, and K is the influence coefficient matrix, which comes from the Green function of Boussinesq’s solution for a point force (). Note that the thickness of the reflective layer should be accurately determined. The small thickness parameter z0 → 0 will induce R → 0 at the original point, which induces a numerically large value K. Thus, the calculated traction force will be underestimated if the thickness of the reflective layer is neglected.

The problem to calculate p from the measured displacement u can be achieved by an inversion process in Fourier space (, ). It is a computational efficient method but matrix K can be ill-conditioned and direct inversion not always work well. In particular, the matrix with high wave number in the Fourier space usually has high condition number, which means the high spatial frequency components of displacement error can be strongly amplified during the inversion process. Therefore, the inversion results notably depend on the appropriate low-pass filter.

Here, we used an iteration method to find the solution p from K * p = u. This method is inspired by the similar approach used for solving the problem of rough contact (, ). We extend it to calculate the deconvolution. To find the solution p of K * p = u, we define a cost function of the force field W(p), and the problem is transformed to minimize itThe discrete convolution of matrix is defined as (x * y) = xy. The symbol “⨀” in this work is defined as the element sum of matrix after element-wise multiplication, which takes the form of x⨀y = xy. It can be proved that if W(p) get minimum at p*, then . Thus, we can get the solution p* of K * p = u. The deconvolution process can be achieved by a gradient decent iteration procedure (more details in section S1). Note that the Fourier transform and inverse Fourier transform can be used to speed up forward convolution calculation in each iteration.

When the deformation or the characteristic traverse contact scale is comparable to the substrate thickness, the effect of substrate thickness requires to be considered, and the SI assumption does not hold. A mechanical model that considering the limited thickness (LT) of substrate was proposed by Álamo et al. () and extended by Xu et al. (). This theoretical model can be explicitly calculated by the matrix product in the Fourier space of the xy plane, where the and are the force and displacement field in Fourier space, h0 is the height of the feature layer, and h is the thickness of the substrate. The matrix Q represents the mechanical response of the substrate with certain thickness [more details can be found in () and ()]. The comparison of the numerical methods to calculate traction stress is discussed in section S2. On the basis of different applicability and advantages, the mechanical model considering the LT is used in the snail crawling force measurement, and the iteration method for the SI mechanical model is used in other experiments.

Validation of the resolution and accuracy

For the demo structure in this study, the spatial resolution can be changed from 100 μm over a 60 mm–by–50 mm field of view to 10 μm over 6 mm–by–5 mm field of view. Meanwhile, the time resolution can be about 11 ms, determined by camera frame rate. They can be readily extended by improving the frame rate and resolution of the camera, as well as providing high-quality feature layer (optimal size of the speckle pattern is 3 to 5 pixels). The stress resolution is determined by the effective elastic modulus of the elastomer, spatial resolution, and the size of contact region. We proposed a basic formula δp ∙ λ~E ∙ δz to describe it, where p is the traction stress, λ is the characteristic contact length, E is the effective elastic modulus of the elastomer, and z is the deformation (more details in section S3). For a determined spatial resolution (δz), using a soft elastomer can improve the stress resolution. On the basis of the formula, we can infer that, for the purpose of characterizing contact fingerprint (typically p ~ 10 kPa, λ ~ 100 μm, and δz ~ 20 μm), the appropriate elastic modulus should be at the level of about 50 kPa. On the basis of the method, the 3D fingerprint was successfully displaced when the finger is pressing on the surface (shown in Fig. 2A). The spatial accuracy can be evaluated by a pin pressing test, in which a pin was pressed on the surface and then lifted and moved with a step of 10 μm. The displacement of the contacting center between steps can be obtained, which was 9.96 ± 2.21 μm (shown in Fig. 2B). The force accuracy was evaluated by a ball pressing process based on Hertz contact theory, because the available “ground truth value” is only the total force obtained by commercial force sensor. We use both our iteration method based on SI mechanical model and the LT mechanical model to calculate the traction stress (shown in Fig. 2C). Both deformation and normal stress generally matched well with the Hertz model prediction when deformation is small. When further increased the indentation (e.g., greater than 2 mm, 10% of the total thickness), the effect of substrate cannot be neglected, and Hertz contact model overestimate the deformation and underestimate the pressure stress. The calculated stress using the SI model also underestimated the maximum stress compared with the LT model. Therefore, LT model is supposed to be adapted when the deformation is not much smaller than the total thickness.

Fig. 2.

Verification of the resolution and accuracy.

(A) Reconstruction of the contact geometry of a human finger. (B) Evaluation of the spatial resolution and accuracy by a pin pressing test with a stepped lateral movement of 10 μm. The inset showed the measured displacement between steps (9.96 μm), and the dashed region indicated uncertainty (2.21 μm) evaluated from the SD. (C) The measured depth and maximum pressure when a ball with a diameter of 12.7 mm presses on the surface (with alcohol to eliminate adhesion). The traction stress was calculated on the basis of SI mechanical model and LT mechanical model, which was used to compare with the Hertz contact theory.

Imaging of the adhesive contact

Because of the ubiquity of intermolecular forces, the dynamic contact process usually involves adhesion (). The adhesive contact between objects is basic physical phenomena and has great applications in various fields such as surgical treatment () and gecko-inspired robotics (). Great progress has been made to model the adhesive stress distribution for different contact problems (), but very few methods are available to directly measure it. The proposed method can be effective to characterize the adhesive stress distribution and evolution. First, the validation was evaluated by the theoretical results describing a ball pulling-off process developed by the Johnson-Kendall-Roberts (JKR) model (). The normal stresses were compressive over the entire contact zone as the sphere was pressed on the surface (Hertz state), and tensile stresses were observed near the edge of the contact zone when close to the pulled-off point (JKR state). Note that the uploading process need first underwent a transition from the Hertz contact state to the JKR contact state. The predicted curve based on the indentation d according to Hertz contact model or JKR model is shown in Fig. 3A (details in section S4), which showed a good qualitative agreement. Some biases occurred at the edge of the contact region for the JKR contact state probably due to three reasons. First, the sharp deformation at the edge was filtered by the optical system considering the limited resolution, and thus, the deconvolution procedure was also not very sensitive to sharp boundary. Second, the polymer elastomer is naturally viscoelastic, which is somewhat biased from the theoretically predicted pure elastic response. Third, it somewhat reflected the controversial point of the JKR model, which theoretically predicted an infinite stress at the edge. In practice, this stress singularity is suppressed by the limited yield strength of material or a finite adhesion stress at a cohesive zone instead of infinite stress at a line ().

Fig. 3.

The normal stress for adhesive contact.

(A) The force-time curve and stress distribution when a Si3N4 ball with a diameter of 12.7 mm was pressed and then pulled off from the sensor surface. The effective work of adhesion W is about 2 mJ/m2, determined by pull-off force based on JKR model. The Hertz or JKR fitting was predicted on the basis of the indentation. (B) The force-time curve during a micropatterned surface contacting and separating from the substrate surface. The stress distribution was a mirror display for better visualization. For both tests, the measured force was obtained by integrating the normal stress. The ground truth force was recorded by a commercial force sensor.

Micropatterned surfaces are believed to enhance adhesion and therefore draw continuous attention (). Here, a silicone surface modified with a gecko-inspired microcylinder array was used to perform adhesion testing. During the contacting process, the microarrays were contacted first, and the pressure stress increased gradually. When it was pressed deeper, even the substrate also contacted with the surface. During the separation process, the normal stress decreased with two apparent steps in the force curve. The first step occurred when the silicon substrate detached with adhesion stress occurring at the edge region, while the second step occurred when the microcylinder array totally detached (shown in Fig. 3B and movie S1). Imaging the stress evolution during contact and detach process would provide important implication for microarray designing, reliability monitoring, and adhesion controlling in application. Moreover, it showed the extraordinary ability of the method to discriminate the multipoint contact state even under overlapped contact forces.

Imaging of rolling friction contact

The ability to discriminate both normal and tangential stress can provide detailed information of interfacial friction (). Rolling is the most common and effective way to reduce friction. Up to now, several qualitative and quantitative rolling friction theoretical models have been proposed (–). However, traditional experimental methods can only provide indirect verification of the total friction force, even with the help of a high-precision force sensor. The proposed method can provide unique information for understanding rolling friction.

The experiment system is shown in Fig. 4A. A ball was placed between two silicone plates. The upper silicone plate was fixed on a commercial force sensor to record the ground truth force. The lower one was designed on the basis of the proposed method to record the interfacial stress. It was dragged by a linear motor stage to drive the ball rolling. Ethanol was used as a lubricant to reduce the interfacial adhesion to compare with dry condition (the slip-roll ratio is measured to be small in section S5). During the loading process, the normal stress is a bell-like symmetric shape and is consistent with the Hertz contact theory for both dry and lubricated conditions. The lateral stress (along rolling direction) and transversal stress (perpendicular to the rolling direction) presented an antisymmetric distribution, and the summation was nearly zero. After the ball rolled, the rolling friction of lubricated condition is slightly smaller than that of dry condition (shown in Fig. 4B). This slight difference can be readily understood in the stress distribution (shown in Fig. 4C). For lubricated condition, the resistant lateral stress in the front-side region became greater than the thrust lateral stress in the rear side, and the difference determined the lateral friction force. While for the dry condition, an adhesion stress region occurred at the rear edge and induced additional stress for lateral friction (the movie is available in movie S2).

Fig. 4.

The contact stress for rolling friction.

(A) Schematic illustration of the rolling friction test. A Teflon ball was placed between the polydimethylsiloxane (PDMS) plates. The load is 0.1 N and the driving velocity is 2 mm/s, corresponding to 1 mm/s of the rolling velocity (more details in section S5). (B) The measured normal force curve and lateral force curve with and without ethanol lubrication. (C) The normal stress and lateral stress (along the rolling direction) during dry or lubricated rolling for the moment indicated by the dashed line in (B).

From the view of energy dissipation, viscoelastic hysteresis and adhesion hysteresis dominated rolling friction dissipation (, ). These two components were clearly displayed at different contact regions in our experiment. Note that, although the proposed method is based on the theory of elasticity, the measured forces by integrating contact stress showed an unexpected agreement with the ground truth force. It implies that the viscous term in viscoelastic material is the origin the asymmetric deformation during rolling. Once the asymmetric deformation was formed, the elastic force dominated the resistance, which is consistent with assumption of the stress distribution in some rolling friction models (). Systematic experiments based on the proposed contact stress–measuring method in future work can provide direct and important experimental results for understanding rolling friction and the development of theoretical models. Moreover, because the process of charge transfer is closely related to the state of frictional contact at the interface, the proposed method is expected to have potential applications in contact electrification or triboelectricity.

Imaging of crawling stress of snails

Snails are endowed with the unique crawling ability with only its single-belly foot. The train of muscle deformation pulse traveling along the foot’s ventral surface known as the “pedal waves” importantly participates in crawling (). The most notable and widely accepted interpretation for the basic crawling mechanism is “adhesive locomotion,” which claimed that the mucus under the foot can undergo a yield-heal cycle during shearing and, thus, it can reduce the forward resistance without lifting the foot (, ). However, one subject still left to be clarified is how the normal force participates in the adhesive locomotion, which is key for snails crawling vertically or upside down. Obtaining the 3D traction stress between the snail’s foot and substrate can provide key information for understanding the locomotion mechanism (, ) and may inspire new paradigms of soft robots (–).

A snail (Achatina fulica var.) was allowed to freely crawl on the horizontally placed surface. The normal and tangential force distributions were recorded in situ and in real time (shown in Fig. 5 and movie S3). Note that the traverse contact scale of the snail is comparable to the substrate thickness; the LT mechanical model was adopted to calculate the traction force. The wave, interwave, and rim regions on the snail foot surface can be distinguished from force distribution. In general, the propulsive tangential force generated at the interwave regions of a snail was offset by the resistant tangential force at the rim and waves. Positive (downward) normal force was mainly generated under the rim and anterior part of interwave region, whereas the negative normal force was generated under the wave region and posterior part of the interwave region. The interplay between normal force and tangential force can provide essential information for understanding the crawling mechanism.

Fig. 5.

Crawling stress of a snail and the mechanical model.

(A) Picture of the snail in the experiment and the typical wave, interwave, and rim regions. (B) The measured normal and lateral stress distribution during the snail crawling. The arrows indicate the moving direction. (C) Schematic of the crawling model for a snail. The soft substrate was deformed upward along with the wave region, which generated negative pressure due to the contraction of the oblique muscles. The wave region propagated forward under the pressure of posterior end of the wave, which can easily peel off the interwave region forward.

The segmented negative (upward) stresses suggested the suction effect of wave region and sealing effect of the rim region, which may come from the contraction of muscle in both lateral direction and vertical direction. Considering the incompressibility of the mucus, as well as the large viscous stress penalty for lifting indicated by Denny (), the thickness of the mucus can be negligibly increased at the concavity region. Instead, the soft substrate would deform upward in unison with the foot (). As long as the concavity shape is formed, it will be pushed forward by the hemocoelic hydrostatic pressure generated behind (). This is indicated by the measured high pressure generated at the anterior edge of the interwave. The concavity wave can propagate forward easily as a “rolling-in” and “peel-off” behavior (). This kind of movement on a soft substrate can also reduce the forward shearing resistance of the wave region.

Moreover, we inferred that, in addition to the viscoelasticity of mucus, the pressure-regulated friction also contributed to the adhesive locomotion mechanism. The flow-like mucus was sealed in the wave region with a negative pressure, whereas the interwave region was counteracted by high-pressure stress and, thus, it was under a boundary lubrication regime separated with solid-like mucus. The friction force generated in the interwave region was determined not only by the shear rate but also by the normal pressure. An indication is that, when a snail was removed from the substrate, if it was unwilling to leave, the snail would increase the suction force to generate large-pressure force on the interwave or rim regions as a resistance (the movie is available in movie S4).

Snails also have the versatility of crawling on horizontal planes, vertical walls, and even ceilings. The collaborative mechanism of adhesion and friction force can be clarified according to the 3D stress distribution. A small snail [Bradybaena (Acusta) ravida] was allowed to crawl freely on surface with different inclinations to record the crawling traction stress (results shown in Fig. 6A). In addition to the pedal stress on the small scale, the total stress was modulated by a normal stress on a large scale. Here, we mainly focused on the large-scale stress distribution. On the horizontal plane, the snail was attached on the substrate similar to a suction cup. That is, pressure stress was in the outer region and negative pressure was in the inner region. On a vertical wall, the anterior part of the snail foot tends to generate a negative pressure, and the posterior part tends to generate a positive pressure so as to balance the gravitational moment. On the ceiling surface, the snail can crawl upside down, and the negative pressure can be generated by any part of the snail foot. Thus, it was supposed that the snail’s foot was composed of suction cups with different scales. The small-scale suction cups constituted the pedal waves to propel the body, whereas the large-scale suction cups can provide segmented sealing and sucking capability so that the body can fully or partially adhere to the wall or ceiling.

Fig. 6.

The crawling ability for a snail on different inclinations.

(A) The 3D contact stress when a snail crawled horizontally, vertically, and upside down. The dashed lines outline the approximate shape of the snail. The open arrows indicate the motion direction. The arrows on the upper right corner indicate the positive direction of traction stress. (B) The sectional normal stress evolution map obtained from the stress distribution with different time series. The dashed arrows indicate the stress evolution for a fixed point on the ground (POG) or a fixed point on the snail foot (POS). The slope of the stripes indicates the velocity of the snail (VOS) or pedal waves (VOW).

A normal stress evolution map can be obtained if we averaged the stress distribution in a line along the moving direction and arranged them along y-axis direction with different time series. Take the vertically crawling for example (shown in Fig. 6B), the stripes contained the motion information of the snail. Because the anterior part was kept positive on the substrate, the slope of this portion represented the moving velocity of the snail body (VOS = 0.53 mm/s). Whereas in the center region of the snail foot, the pedal wave propagated much faster and showed a higher slope (VOW = 3.93 mm/s). Moreover, the multiscale stress features can be readily recognized in the stress map. If we focused on the fixed point on the ground (POG; vertical section line in the map), then the periodical stress was modulated by a negative pressure and then positive pressure. However, if we focused on the fixed point on the snail foot (POS; inclined section line parallel to the body velocity line), then the periodical stress kept steady with time. We inferred that the small-scale crawling locomotion and large-scale sucking adhesion is realized by independent body structures of a snail. The findings can provide inspirations for designing new paradigms of soft crawling robots.

DISCUSSION

This work proposed a novel method to characterize the dynamic 3D distribution of traction stress with high spatial and temporal resolution, which can overcome the shortages of the existing methods listed in table S1. This was achieved by combining the binocular stereo vision techniques and elastic mechanical models. Note that the LT mechanical model is theoretically accurate and computationally efficient when the deformation or characteristic traverse contact scale is comparable to the substrate thickness, whereas the introduced iteration method can substantially improve the numerical stability and reduce the dependence of low-pass filter parameters when dealing with the SI mechanical model.

By means of its advantages in characterizing both normal and tangential stress distribution, the proposed method was demonstrated to be powerful for studying many basic and essential problems in interfacial phenomena. For example, the adhesion stress evolution process of a micropatterned surface was visualized, which can be helpful to the design functional microstructures. It was experimentally displayed the contribution of elastic resistance and adhesion traction to rolling friction at different contact regions. On the crawling mechanism of snail, we disclosed that the distributed sucking and sealing effect of the concavity pedal waves with roll-in and peel-off behaviors propelled a snail crawling, and the multiscale suction and sealing effect was responsible for its versatility ability of crawling in the horizontal, vertical, and upside-down directions. We believed that the proposed dynamic temporal-spatial traction stress imaging method could be beneficial to various fields such as tribology, biology, bionics, and robotics.

MATERIALS AND METHODS

Preparation for the measurement system

The picture of the experimental system was shown in fig. S1. It mainly consists of an optical system and an elastomer system. The core optical device was of a commercial binocular camera (CAM-AR0135, Wuhan LenaCV Corp. on Taobao website, China), which can read RGB (red, green, and blue) image frames up to 90 fps at a resolution of 640 pixels by 480 pixels with a USB (universal serial bus) 3.0 connection to a computer. The two cameras were fixed at an angle bias of about 15° from the parallel direction. The elastomer system was designed with several layers. The bulk transparent block was obtained by pouring polydimethylsiloxane (PDMS) silicone gel (Dow Corning 184; Young’s modulus can be adjusted between 0.6 and 2.0 MPa with different cross-linker agent ratios; 9400B, Shenzhen Hongyejie Tech. Co. Ltd.; Young’s modulus is about 20 to 60 kPa). We mixed them together to obtain the desired Young’s modulus into a 3D printed box mold. Two parallel glass plates were inserted in the mold to guarantee the desired thickness of the elastomer, as well as the surface parallelism and smoothness. The thickness of elastomer is controlled ranging from 10 to 25 mm, which determined the maximum deformation range. After the silicone was cured, random speckles were transferred to the surface by two methods. Small speckles (~10 μm) are obtained by an airbrush, and medium speckles (~200 μm) are obtained by a preprinted transfer sticker. A thin reflective layer (thickness of about 80 to 120 μm), which is fabricated by spin-coating the silicone painting and then cured, was covered on the top surface. The silicone painting was prepared by adding titanium dioxide pigment (mass ratio of about 5%) into the same kind of silicone gel.

Procedures to obtain the deformation field

The obtained speckle images of the deformed surface were processed first by the DIC algorithm. In this work, the 3D-DIC technology was achieved mainly based on MultiDIC (), an open-source software tools in MATLAB platform (The MathWorks Inc., USA), as well as the embedded 2D-DIC open-source software Ncorr (). The algorithm details can be found in () and (). Briefly, the procedures to perform the 3D-DIC measurement were as follows: First, two cameras view a region of interest from different positions and capture a sequence of images. Next, 2D-DIC was performed between a given set of points and matching these points throughout the sequence of all stereo images. Then, the correlated sets of image points were used to reconstruct 3D surface shape and 3D displacement field on the basis of precalibrated stereo parameters. The deformation field can be obtained from the displacement field by subtracting the rigid body displacement, followed by some preprocessing procedures including coordinate correction and image filtering. Then, the deformation fields can be used to calculate the contact stress from the elasticity theory. The details to calculate the stress field can be found in section S1.

Preparation for the snail’s experiment

The small snails [B. (A.) ravida] were obtained from the garden in Tsinghua University. Their weights ranged from 3 to 8 g, and their body lengths ranged from 25 to 40 mm. The big white jade snails (A. fulica var.) were obtained from a pet store on Taobao website. Their weights ranged from 20 to 30 g, and their body lengths ranged from 40 to 60 mm. The snails were fed daily with fresh water and vegetables at the proper temperature and humidity. Considering the small contact pressure, the low elastic modulus of the elastomer bulk should be used to improve sensitivity.