Novel Contact Modeling for High Aspect Ratio Soft Robots

Contact modeling between a soft robot and its environment is challenging due to soft robots' compliance and the difficulty of embedding sensors. Current modeling methods are computationally expensive and require highly accurate material characterization to produce useful results. In this article, we present a contact model that utilizes linear complementarity and Hencky bar-chain methods, and requires only static images to efficiently predict the interaction between actuator and environment. These methods have yet to be introduced to the soft robotics community for modeling robots that deform due to eigenstrains or strains not caused by external forces. We validated our model using a custom experimental setup and computer vision algorithm on 3-mm OD, 90-mm long, tube-like actuators. Our results indicated a 1.06% difference in shape between model and experiment, with computation times in 10 s of ms—three to four orders of magnitude faster than nonlinear gradient descent. Additionally, the error in interaction forces between the model and experiment decreased as pressure increased, with an average error magnitude of 45% and 21% for pressures at the low and high ends of the tested range, respectively.

the benefits of having an inherently safe robot with variable stiffness come with tradeoffs. The pliancy of soft robots makes modeling contact (i.e., shape change and interaction forces) with other objects nontrivial. Yet, such information remains critically important in situations where such knowledge is vital to safety.
One example of where the knowledge of these two aspects, shape change and force, is imperative is in minimally invasive surgery. Surgeons must navigate their robotic tools while interacting with delicate anatomy. To ensure patient safety, the contact forces between the robot and its environment must remain below a safety threshold while also being high enough to perform effectively in situations like clearing plaque, anchoring within arteries, or manipulating tissue. In many cases involving small anatomy, such as neurosurgery, researchers seek to reduce the size of their robot to millimeter scale [1]. The reduction in size drastically reduces torque and force transmission along the length of the robot. To overcome these challenges, a contact model that accurately predicts the constrained shape and reaction forces generated over the contact area between a soft robot and its environment is necessary.
The core contribution of this work is a theoretical contact model and experimental validation of the shapes and forces generated by soft actuators when constrained by a rigid environment. The model applies Hencky bar-chain and linear complementarity methods to soft robots to model both self-deformation and contact modeling. Although these methods have been used for other applications outside of robotics [2], [3], the combination of these methods have yet to be applied to soft actuators that deform due to eigenstrains; that is, deformation caused not from external forces but from internal actuation forces (e.g., hydraulic pressure) and the structural design of the actuator. Moreover, modeling using the proposed work can be completed using only static images, thus eliminating the need for embedded sensors, real-time measurements, and detailed knowledge of the actuator's material properties. Ultimately, the proposed work benefits the soft robotics community by introducing a contact model that is generalizable to any class of soft actuators that deform due to eigenstrains and the model does so using computationally tractable algorithms that can be computed in near real time.
Other contributions of this work include a comparison between the linear complementarity approach and nonlinear gradient descent, as well as a computer vision algorithm to determine the centerline of soft actuators, both unconstrained and constrained, as a function of input pressure. B. Related Work 1) Hencky Bar-Chain: The Hencky bar-chain approach to model the bending and buckling of slender structures was first introduced in 1920 by Prof. H. Hencky [4]. This method discretizes structures into a series of rigid bars and torsional springs. Since 1920, many researchers, particularly in the fields of Civil and Structural Engineering, have expanded on Hencky's work, which is summarized comprehensively by Wang et al. [2].
Some examples of the Hencky bar-chain applied to engineering problems include the buckling of columns and frames subject to external axial loads [5], the bending, buckling, and vibration response of elastica [3], the dynamics of uniformly loaded elastica [6], the buckling and vibration in Euler beams with various end conditions [7], and the dynamic response of a beam colliding with rigid obstacles using the Hencky bar chain [8]. However, the listed examples focus on the deformation of beams and elastica under external loads rather than self-deformation as discussed in this article.
The Hencky bar-chain approach is reminiscent of the pseudorigid body (PRB) model, which is often used to model compliant mechanisms, in that it approximates structures using both rigid components and torsional springs [9]. However, the Hencky bar chain does not rely on equivalent parameters to the PRB, nor does this work focus on a mechanism comprised of multiple actuators. Furthermore, there is no restriction on magnitude of the displacement with the Hencky bar-chain method, which has been used to solve a variety of elastica problems [2]. Leveraging these benefits, the work in this article can be applied to a wider variety of actuators beyond those actuated using fluid power.
2) Linear Complementarity Method: The linear complementarity method was first proposed in 1966 by R. Cottle as a mathematical optimization problem [10]. Given a matrix K and a vector z, the solution of the problem seeks vectors R and g that satisfy the following constraints: This method is used in modeling the contact and spacing between an object and a constraint given a complete description of the system.
Specific to soft robotics, several groups have formulated their problems as linear complementarity problems (LCPs) [11], [12], [13]. Al-Fahed et al. [11] focused on exploring the difference between two finger grippers, "hard" and "soft," and found that soft grippers have the added advantage of constraining the rotation of objects. Ciocarlie et al. [12] use known expressions for contact between elastic bodies to solve an LCP to model soft fingers.
In medicine, Katz and Givli [14] examined the post-buckling behavior of a beam constrained by springy walls to model how guidewires deform when inserted into an occluded artery. However, the work differs in that the deformation of the beam was dictated by an external axial force and utilizes energy methods to find a solution.
Many have modeled the contact mechanics for soft or continuum manipulators [15], [16], [17]. However, these models also require an external, axial force to describe contact and do not use the same approach as proposed in this article. Yip and Camarillo [18] used a unique, model-less approach to control tendon-driven manipulators in constrained environments, but primarily focused on closed-loop control rather than the quasistatic estimation of soft robot poses and forces in this work. The work of Shiva et al. [19] is directly relevant to the proposed work in that it produced a quasi-static model for position and contact force estimation of soft robot appendages in sensor-deprived environments, but differs from this work in the methods used.
Other groups have utilized the deformation of soft robots to locomote in both flat and tube-like environments, including in [20], [21], and [22]. These groups and others, including in [23], have utilized soft robotics' pressure-controlled curvature and flexural rigidity to explore how those variables affect behavior and performance.
Finally, finite-element method (FEM)-based methods have been commonly used to model deformation of soft robots [13], [24] as well as contact problems more generally [25], [26], [27]. While FEM methods may provide rigorous modeling, the drawbacks of FEM approaches include the large computational cost and requirement for detailed material properties for each actuator component [24], [28], [29]. Pozzi et al. recognized this tradeoff between computational expense and accuracy, and proposed a variation on traditional FEM using kinematic chains to reduce the computation time, but the method has drawbacks in not being able to predict nonlinear movements such as buckling [30].

A. Problem Description
To accurately describe the contact forces and shape of a soft actuator when constrained by its environment, we must first describe the unconstrained shape of the actuator and its relation to input pressure, P . The actuator used throughout this article is representative of a popular design used in soft robotics [31], and is comprised of three key components, shown in Fig. 1(a): base tube made from elastic rubber, strain-limiting layer, and inextensible fibers. The angle at which the fibers are wrapped generates extension when the actuator is pressurized, which makes the actuator a fiber reinforced elastomeric enclosure (FREE) [32]. Conversely, the strain-limiting layer prevents extension along the side of the actuator on which the layer is placed. This causes each strain-limited section of the actuator to bend. This tendency for self-deformation can be modeled as an eigenstrain ("eigen" in German means "inherent" or "self"), meaning the deformation is not caused by external mechanical stresses.
Multiple strain-limiting layers placed in an alternating pattern along the length of a base tube form a planar anchoring actuator as shown in Fig. 1(b), which can be used to anchor within or locomote through tube-like structures [22]. Each segment section within the anchoring actuator can be described by its overall length, which is the same length as the strain-limiting layer, as well as its bend angle, θ, and bend radius, r, shown in Fig. 1. The bend angle, θ, in each section of length can be expressed as where θ increases with input pressure P [see Fig. 1(c)]. Our model will rely solely on bend angle θ, and we will later derive an experimental relationship between θ and P that will change based on material properties of the actuator. When the actuator is constrained by its environment (e.g., anchored into an artery or pipe), the actuator assumes a new shape and reaction forces are generated along the contact area. Both of which are dictated by actuation pressure and the height of the environment, shown in Fig. 1(d).

B. Notation for Unconstrained Actuator
Throughout this article, scalars are represented using lowercase letters not in bold typeface (e.g., d, μ, and b 1 ). Vectors are lowercase, bold typeface letters (e.g., α,τ , and c 1 ). Matrices are uppercase, bold typeface letters (e.g., A, B, and L). The exception to this is the reaction force vector,R, which is one of the variables for which we are solving. A bar above the variable  Fig. 1(c). (b) Discretized description of soft actuator of N links and N + 1 nodes, where N = 6 in this example. Note that the furthest right node is labeled as n for notation simplicity, where n = N + 1. Each link is subject to self-weight, μ. In the unconstrained configuration, the orientation of the ith link relative to the positive x-axis is defined as α * i , φ * i describes the orientation of the ith link relative to the i-1th link, and u * i describes the displacement of the ith node in the y-direction. In the unconstrained configuration, the torque at the ith node, τ * i , is zero. When the actuator is constrained, each node undergoes an angular displacement and a nonzero torque, τ i , is developed.
(e.g.,φ andū) represents that the vector has been truncated to only consider the variables within the vector definition (e.g., [u 2 , . . . , u N ]) instead of variables 1 through N . We use the following notation to describe an unconstrained actuator: α * is an N × 1 vector describing the ordinary angles between each link and the horizontal,φ * is a (N − 1) × 1 vector describing the angles between each pair of links (joint angles), u * is a (N − 1) × 1 vector describing the vertical displacement of each node, andτ * is a (N − 1) × 1 vector containing all zeros for the unconstrained actuator. This is because there is no unconstrained torque due to the "eigen" curvature of the actuator, where the unconstrained curvature is induced by fluid pressure and not external forces. Each variable is shown in Fig. 2 and the star * notation denotes the unconstrained "eigenstrain" configuration.
Our unconstrained actuator description is based on an Nlink discretization of the actuator bending section described by Fig. 2 based on the Hencky bar-chain method. Each of the N + 1 nodes can be thought of as an elastic rotational spring capable of generating torque when the rigid links are displaced. The actuator end conditions can be described as pinned-pinned, or simply supported. As such, the first and last nodes do not undergo any vertical displacement, nor do torques develop at the first and last nodes, because of these boundary conditions. Fig. 3. Iterative process for determining the final solution for reaction forces and actuator configuration. A refers to the unconstrained configuration ( * notation) found using the Hencky bar-chain method, B represents the previous, known configuration used as a "trial" value in the iterative process, C represents the constrained configuration after solving each subsequent iteration of the LCP (∧ notation), D represents a conditional statement used to determine whether the maximum change in link orientation (max(|α − α|)) is less than a threshold, E refers to the intermediate configuration, and F is the final configuration.
We use this unconstrained description in two ways. The first is to create a baseline, eigenstrain configuration on which we impose the environmental constraint. This is further outlined in the next section. The second is to later derive a relationship between actuation pressure and the resulting bend angle.
The following are assumed.
1) The links used in the Hencky bar chain are rigid.
2) The torsional springs connecting the links are frictionless.
3) Contact between the actuator and its environmental constraint is frictionless. 4) No torques have developed at the nodes in the initial, unconstrained configuration. 5) The actuator is uniform, meaning: a) the elastic modulus, E, as well as stiffness, EI, are known and constant across the pressure range of interest; b) the weight per unit length, w, of the actuator is constant.

A. Problem Definition
As described previously, the major contribution of this article is a description of the actuator shape and the forces developed along the contact area when the actuator is constrained by its environment. This is achieved by discretizing the actuator using the Hencky bar-chain approximation and formulating the constraint problem as an LCP. More specifically, we linearize the equations describing the actuator then solve a sequence of linear problems corresponding to a changing configuration of the actuator to ultimately reach a final nonlinear solution. The process is shown in Fig. 3 and the notation is outlined in Section III-B.
Due to the assumption of frictionless contact between the actuator and its environment, as well as the pinned-pinned end conditions of the actuator, we are only concerned with forces in the y-direction.

B. Notation for Constrained Actuator
Throughout this article, the hat notation is used to communicate when a variable represents the unknown configuration of the constrained actuator (e.g.,α,φ, andū) and a nonhat i represents the y-direction displacement of each node and g * i represents the y-distance between each node and the constraint. (b) Constrained discretized description of a soft actuator, where theĥ at notation specifies the constrained orientation.
to represent the previous, known configuration of the actuator (e.g., α,φ, andū). In addition to the hat and bar notation, as well as the variables introduced in Section II-B, we use the following variables to describe an actuator constrained by its environment:ḡ is a (N − 1) × 1 vector describing the vertical gap between each node and the environment,h is a (N − 1) × 1 vector representing the height of the environment, andR is a (N − 1) × 1 vector describing the reaction force at each node.
These variables, shown in Fig. 4, are used in the LCP formulation to solve for the reaction forces,R, and the constrained configuration using the gap between each node and the environment,ḡ.
In summary, the inputs to the model are as follows. 1) Geometry: Overall length, , bend angle, θ, and constraint height, h. 2) Actuator material properties: Weight per unit length, w, and bending stiffness, EI, which may be empirically determined [33]. 3) Discretization resolution: Number of links, N . The outputs of the model are as follows: reaction force at each node,R, and y-direction gap between each node and the constraint,ḡ.

C. Linear Complementarity Problem Formulation
Formulating the constrained actuator problem as an LCP allows us to not only determine the reaction forces developed along the length of the actuator created through contact with the environment, but also to understand the constrained configuration (e.g., position of each node and angle of each link). Understanding both of these variables is critical in delicate environments and situations where we must know how the actuator shape change and interaction forces change to ensure adequate contact between the robot and its environment (e.g., traction for anchoring, locomotion, and tissue manipulation) without exceeding force limitations.
We must first determine expressions for the support reactions q 1 and q n as functions of the external reaction forces,R. Summing the moments about node n generated from external forces, we get where d is defined as The weight of each link, μ, acts at the center of mass and is defined as We can express b 1 as We can define c 1 as an (N − 1) × 1 vector Following a similar process, we obtain the reaction force at node n by summing the moments about node 1 generated by external forces. We obtain The expression for b n is and c n is a (N − 1) × 1 vector defined as We can calculate the torque acting on each node by summing all of the moments generated by external forces, resulting in τ = a 2 (γ + AR) (11) where a represents the length of each link, defined as and γ is a (N − 1) × 1 vector defined as (16) shown at the bottom of this page. We can define a second expression for torque by relating torque and the change in joint angle from the unconstrained configuration,φ * , using the constitutive equation for a torsional spring.τ where EI represents the bending stiffness of the actuator. Utilizing kinematics, we can express the relationship between φ andα asφ = Lα (18) where L is an To relate small y-direction displacements of each node,û, to the angular deflections of each link,α, we must first linearize Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.
about the previous value for α using The linearization results in the following expression forα: and the displacement of each node is a (N − 1) × 1 vector and Solving forR, we must first combine (18) and (21) to get an expression forφ asφ Combining (17) and (25), we can express the torque at each node as a function of the y-direction displacement of each node.
Comparing (11) and (26), we can solve forR as (27) The relationship between the y-direction displacement of each node and the y-direction gap between each node and the environmental constraint can be expressed as (see Fig. 4) Combining (27) and (28) and solving for the reaction forces at each node,R, we getR where K represents the system's stiffness matrix, expressed as and We can then use the expressions for K and z to solve the LCP (e.g., LCPSolve in MATLAB, [34]). To find a nonlinear solution to the contact problem, we solve the LCP iteratively until The threshold of 0.001 radians was chosen to ensure the LCP has converged to a solution with minimal change in link orientation. The final, nonlinear solution provides the reaction forces at each node,R, as well as the gap in the y-direction,ḡ, between each node and the environment constraint. We can useḡ to deduce the actuator's constrained shape.
With the proposed problem formulation established, we now discuss how the formulation was implemented to perform analysis on the relationships between the bend angle, length, and material properties of an actuator and the shapes the actuator obtains and the reaction forces developed when constrained by a rigid environment.

D. Determination of Discretization Resolution
The number of links used for discretizing the soft actuator was determined using a thresholding method. We solved the LCP to find the total reaction force produced for actuators with unconstrained bend angles of 10 • through 180 • (10 • intervals), constraint height of 2 mm through 7 mm (1-mm intervals), and a constant bending stiffness and actuator length of EI = 32.01 N · mm 2 and = 30 mm. Note that these values were chosen to align with the experimental work presented in Section IV, where EI was found experimentally [33]. For each combination of bend angles and constraint heights, we increased the number of links in the discretization by 10. A sufficient number of links was determined when the percent change in total reaction force was below a 0.01% threshold.
A conservative estimate was used when determining the necessary number of links, meaning we used the largest number of links until every percent change in reaction force throughout the entire space was below the threshold. These results are shown for h = 3 mm in Fig. 5.
Further analyses throughout Section III utilize discretizations of 140 links as 140 links yielded sufficient discretization for the range of bend angles listed previously.

E. Theoretical Shape of a Constrained Actuator
The relationship between the constrained shape of the actuator and its unconstrained bend angle, θ (e.g., due to an increase in actuation pressure) can be seen in Fig. 6. As the bend angle increases, the actuator comes into contact with the environmental constraint, which generates a reaction force. As the bend angle continues to increase, the actuator's shape flattens and the contact area increases. Fig. 7 shows the relationship between the actuator's shape and the height of the environmental constraint. To effectively explore this relationship, the bend angle was kept constant.
1) Observations: Theoretical Shape of a Constrained Actuator: We observed that the actuator has a tendency to flatten when in contact with a rigid constraint. The amount of flattening or the amount of actuator in contact with the constraint, increased both with a decrease in constraint height as well as an increase in bend angle. The results also showed that the constrained actuator span (x-direction distance between the ends of the actuator) increased in comparison to the unconstrained actuator span, where a decrease in constraint height corresponded to a larger increase in span. This was because the overall actuator length stayed constant, but vertical deflection was limited by the constraint.
The model provides a better understanding of the relationships between actuator shape, bend angle, bending stiffness, and constraint height to use as both a design tool and predictive model of the actuator contact behavior. This understanding of actuator shape could be used to predict surface traction and contact, estimate locomotion performance by understanding how the actuator's span changes as it is pressurized, and avoid critical structures by understanding the actuator's configuration.

F. Theoretical Reaction Force of a Constrained Actuator
The relationship between reaction force and bend angle was determined by keeping constant the actuator length ( = 30 mm), number of links (N = 140), and bending stiffness (EI = 32.01 Nmm 2 ) while increasing the bend angle from 0 • to 180 • . The results are shown in Fig. 8(a). The regions where the total reaction force vanishes represent the bend angles for which the actuator was not in contact with the environment.
The interaction between reaction force and bending stiffness, EI, was explored by keeping the actuator length ( = 30 mm), number of links (N = 140), and bend angle (θ = 150 • ) constant while increasing the bending stiffness from 10 to 100 N·mm. This was repeated for a range of constraint heights (2-7 mm). The results are shown in Fig. 8 Similar trends were seen after examining how the reaction force versus bending stiffness relationship changes as the constraint height stays constant but the bend angle varies. This was determined by keeping constant the actuator length ( = 30 mm), number of links (N = 140), and constraint height (h = 3 mm) while increasing the bending stiffness from 10 to 100 Nmm. This was repeated for the full range of bend angles (10 • to 180 • ) to ensure contact between the theoretical soft actuator and the environment. The results are shown in Fig. 8(c).

1) Observations: Theoretical Reaction Force of a Constrained Actuator:
The results showed that the largest reaction forces consistently occurred for the smallest constraint heights. As expected, this implies that the more deformation required to reach the constraint, the more force is lost to deformation. We also saw a consistent increase in the total reaction force as both the bend angle and bending stiffness were increased. This implies that stiffer actuators may exert more force when in contact with a constraint. However, in the context of soft robotics, increasing the stiffness of actuators comes with a potential tradeoff of reduced deflection or manipulability (i.e., bend angle in this context). Thus, our proposed model could be used as a design tool to quantitatively model the tradeoff between bending stiffness, force exertion, and kinematic manipulability to determine the most suitable; combination for a given application.
Overall, the model provides a better understanding of the relationships between reaction forces, bend angle, constraint height, and bending stiffness, which could be valuable while actuating in delicate environments where forces must be predicted and controlled. Knowledge of actuator forces could also be used to model robot performance in manipulation or locomotion where interaction with other objects or constraints is critical.

G. Validation of the Linear Complementarity Approach
A cantilever beam was used to validate the LCP approach and explore the accuracy of the approach outlined in Sections III-A and III-C. The cantilever beam was subject to a uniformly distributed load, with the beam a distance h away from an environmental constraint, as shown in Fig. 9. For simplicity, the weight of the beam was ignored.
A similar process to that described in Section III-C was followed and the resulting deflection was compared to the theoretical beam deflection equation expressed as To compare reaction forces, the theoretical reaction force, R t , was determined by first determining the distributed load required, F c , to deflect the beam such that its deflection at the free end is equivalent to h, where The residual distributed load, F r , represents the amount of the load that is in excess of the F c , or the load contributing to the reaction force once the beam is in contact with the constraint, where The deflection at the free end caused by F r can be expressed as . Using Δy r , we can determine the point load, R t , at the free end required to achieve the residual deflection where R t also represents the theoretical reaction force. The results are shown for a cantilever beam of length = 30 mm, EI = 32.01 Nmm 2 , and F = 0.001 N/mm that was unconstrained [see Fig. 10(a)] as well as constrained by its environment [see Fig. 10(b)] and discretized into N = 140 links. The reaction force determined using LCP and a constraint height of h = 1 mm was found to be 7.9 mN and the theoretical reaction force found using beam deflection equations was 7.7 mN, resulting in a 2.60% difference.

1) Observations: Validation of Linear Complementarity Approach:
The comparison between the proposed LCP formulation and the classic beam deflection equation showed that the LCP model was able to predict interaction forces with less than 3% difference to the linear beam theory solution and a normalized difference in deflection of less than 8 × 10 −3 . We observed that the normalized difference in deflection (i.e., Euclidean distance between the LCP and theoretical model scaled by beam length) increased linearly over the length of the beam. A portion of this difference could be attributed to the limitations of the theoretical beam deflections as being valid only for deflections of less than 10% of the overall beam length. Another possible explanation is that small differences in the orientation of the links relative to one another build up over the length of the beam.
The results also demonstrated that the LCP model maintains the expected constraint. That is, the constrained beam deflection in the y-direction never exceeds the limitation imposed by the constraint. Overall, the cantilever beam comparison between commonly used theoretical equations and our proposed LCP formulation reasonably confirmed that the LCP model behaves as expected.

H. Comparison of Two Contact Model Approaches
With the confirmation of the LCP formulation and acknowledging that an iterative LCP approach is not the only method to solve a nonlinear contact problem, we compared the LCP approach to the same contact problem formulated as a system of nonlinear equations. We solved the nonlinear system of equations using a gradient descent optimization approach, which solves a problem specified by F (x) = 0 for the vector x. Further details of the comparison between the gradient descent optimization and LCP approaches are shown in the Appendix.
The comparison showed that LCP appears superior in terms of both accuracy and computation time, particularly in situations where the height of the constrained shape is much larger than the height of the environmental constraint (i.e., large bend angle and low constraint height). The LCP successfully modeled actuators constrained to low constraint heights of less than 10% of the actuator length, provided more information about the distribution of reaction forces, and was significantly more computationally efficient because it solves a linear problem iteratively to find a nonlinear solution rather than optimizing 3N + 2 equations. The comparison showed that the nonlinear gradient descent method requires a normalized CPU time that is of order O(10 3 ) to O(10 4 ) larger than the LCP CPU time, depending on the bend angle (Note: Actual computation time for the LCP approach at 40°was 24.7 ms).

A. Actuator Construction
The actuators used for experimentation in this article were FREEs manufactured using the steps shown in Fig. 11. The base tube shown in Fig. 11(a) was created by injecting liquid, two-part polyurethane (Polytek 74-20) into the custom fixture and mold shown in Fig. 12. The outer diameter of the base tube was dictated by the inner diameter of the borosilicate glass tube placed in the mold. The inner diameter of the base tube was determined by the diameter of the tensioned, stainless steel rod running through the center of the glass tube. The custom fixture and mold enabled the manufacturing of base tubes with selectable material properties and geometries.
Once the base tubes (Inner Diameter, ID = 1.59 mm, Outer Diameter, OD = 3 mm) were removed from the mold, a strainlimiting layer (3-M Durapore tape) was placed in the pattern necessary to generate the desired actuator behavior. The desired behavior to relate our experimental actuator to our theoretical model was a bending actuator, so the strain-limiting layer was placed on one side of the actuator and kept in place using adhesive backing.
The tube and strain-limiting layer were then inserted into a custom CNC lathe, shown in Fig. 13, that wrapped two inextensible fibers (cotton thread) at known angles along the length of the actuator (−78 • , 78 • relative to the longitudinal axis) [22].
The wrapped actuators were then dip-coated in a final, thin layer of polyurethane to prevent the fibers from slipping and make the actuator walls more robust.

B. Free-Fold Bending Stiffness Estimate
Due to the composite nature of our hydraulic actuators, obtaining an accurate prediction of the elastic modulus, E, of the overall structure is difficult. To avoid this, the elastic modulus and area moment of inertia were lumped together into a single bending stiffness term, EI. The bending stiffness was found experimentally using a free-fold test approach, first introduced by Peirce to determine the bending rigidity and bending length for fabric sheets and later adapted for soft robotic actuators by McDonald et al. [33], [35].
Using this method, the average experimental bending stiffness of EI = 32.01 Nmm 2 and weight per unit length of w = 2.45 × 10 −4 N/mm was found for the soft actuators.

1) Custom Test Fixture:
The reaction forces and constrained shape of the soft actuator were tested using a custom test fixture shown in Fig. 14. It is important to note that an experimental actuator comprised of three identical bending sections was used, but the only section of interest for this article was the middle bending section that was in contact with the load cell attachment. The reason for using an actuator comprised of three bending sections rather than one was to minimize the effects that the barbed fittings and input pressure line had on the behavior of the actuator as it was pressurized. The two end bending sections created a symmetric environment for the middle bending section that matched the simply supported end conditions of our model.
2) Data Collection: Data were collected on three actuators. Each actuator was connected to the setup shown in Fig. 14 and pressurized using a syringe connected to the setup via a pressure line. Each actuator was pressurized from 0 to 172.37 kPa (25 psi) in 17.24 kPa (2.5 psi) increments in four configurations: unconstrained and constraint heights of 4.5, 5.5, and 6.5 mm. An image of the actuator was captured at each pressure increment. The pressurization sequence was repeated three times for each actuator for each of the four configurations. Pressurization was done slowly so as to minimize the dynamic phenomena of the actuators.
For the constrained actuators, the constraint heights were determined by changing the separation between the upper and lower constraints shown in Fig. 14(b) H, where the separation Fig. 14. (a) Experimental setup for testing the reaction forces and centerline shapes generated when a pressurized soft actuator was constrained by its environment. Inset: close-up of pin joint used to enforce pinned-pinned end conditions. A barbed fitting was inserted into the actuator to allow pressurization, while two prongs on the fitting were press fit in low-friction bearings to allow the fitting to rotate freely. (b) Front view of setup, where A: Load cell attachment that served as the bottom half of the constraining environment, B: Carriage that was fixed in place to prevent linear translation of the left pin joint, C: Soft actuator, D: Linear guide rail that served as the top half of the constraining environment (PBC Linear, LPM17-0200-1), E: Low-friction carriage that was able to translate linearly to allow axial change in length when the actuator was constrained, F: Bearing holder with adjustable height to ensure the pin joint sat halfway between the environment walls, G: Capped barbed fitting with prongs press fit into the rotational bearing (McMaster, 57155K201), H: Separation between the two plates was set to be equal to double environment height (2h), I: Load cell spacer, and J: Load cell (SparkFun, TAL221, 500 g). (c) Side view of experimental setup, where: K: Nonglare acrylic guidance plates to prevent out-of-plane deflection (McMaster, 1/8" thick, 1178T11). was equal to twice the constraint height (e.g., separation was 11 mm for a constraint height of 5.5 mm). This was done manually using gap gauges and calipers before the experiments. Force data were collected at each pressure increment using a calibrated and tared load cell.
Data were collected at 10 Hz using a microcontroller (Teensy 3.5) connected to a pressure transducer (Honeywell TBP-DANS030PGUCV) and the load cell, each connected to load cell amplifiers (SparkFun, HX711). The camera used for capturing the images (iPhone 11 with Dual 12MP cameras) was placed at a fixed distance from the test setup and at the same height as the actuator.

D. Computer Vision
Once experimental data collection was complete, computer vision was used to determine variables such as the relationship between input pressure and the resulting bend angle of the actuator, length of the bending section, and a quantitative comparison of the similarity between experimental and theoretical centerlines for a given input pressure. Determining the relationship between actuation pressure and bend angle, θ for each actuator was a critical step to compare our experimental force and centerline shape with our theoretical model. This was because the unconstrained bend angle was used as an input to the constrained actuator model, as described in Section II-A. Similarly, determining the length of the experimental bending section was a critical input to the model and allowed us to create a better baseline for comparison.
The desired variables were determined using an interactive program (implemented in Python 3.7.8) to analyze the images captured during experimentation. The analysis consisted of several stages.
1) Preprocessing: All images were rectified before quantitative analysis to account for warping caused by the camera lens (using OpenCV's undistort() function). The images were also leveled by comparing the bearing centers on either side of the test stand. The bearing centers were determined using the circle Hough transform technique as shown in Fig. 15. The difference in the coordinates of each of the bearing centers determined the rotation angle to level the image.  2) Coordinate System Setup: Once the image was leveled, a coordinate system was created using the center of the left bearing as the origin. The pixel-to-mm conversion ratio was calculated using an Augmented Reality University of Cordoba (ArUco) tag of known dimensions (25 mm × 25 mm) fixed to the test stand [36], [37].
3) Segmentation: Segmentation of the image was performed using Hue, Saturation, Value (HSV) color ranging. A preview of the segmentation was displayed and updated as the user manipulated sliders controlling the HSV threshold values. Once a satisfactory segmentation of the image was achieved, the background pixels were set to pure black. The segmentation can be seen in Fig. 16.

4) Centerline Detection:
Using the segmented image, the actuator centerline was detected using an edge-finding technique. Starting in the upper left pixel in the image, the detection was conducted according to Algorithm 1, where each previous pixel refers to the pixel above the current pixel.
The actuator's centerline coordinates were determined by calculating the average of the row numbers of the top and bottom edge pixels for each column in the image. These coordinates were then converted to length units using the pixel-to-mm conversion. Thus, for each position on the centerline, only the y-coordinate had to be calculated, as the x-coordinate was known from the column of pixels currently being analyzed. Example detected edges and resulting centerline are shown in Fig. 16.

5) Critical and Inflection Point Detection:
Using the extracted centerline, the start and end of the middle bending section of the experimental actuator were determined by finding the critical points along the entire centerline. The critical points Algorithm 1. Edge-Finding Algorithm. Fig. 17. Detected centerline, inflection points, and minimum loci for the middle bending section overlaid on the unconstrained actuator. This image is used to illustrate centerline extraction used for the determination of bend angle and inflection points, and is not representative of the typical constrained case.
were defined to be the inflection points (where each bending section starts and stops), as well as the minima and maxima of each bending section. The critical points can be seen in Fig. 17.
6) Circle Fit: Because the middle bending section of the actuator was the section of interest for comparing our experimental results to our model, the relationship between pressure and the bend angle of the middle section was determined using a least-squares circle fit. The circle fit was performed on the centerline between the two inflection points and returned the fit circle's center coordinates and radius [38]. Fig. 18 shows the resulting circle fit for an illustrative sample centerline.

7) Bend Angle and Arc Length:
The bend angle, θ, of the middle section of the actuator was calculated using the center of the circle fit (x c , y c ), as well as the inflection points, (x 1 , y 1 ) and (x 2 , y 2 ), where where The radius of the circle fit, r c , and the bend angle, θ, were then used to calculate the arc length of the centerline points between

8) Centerline Coordinate Comparison:
Quantitative comparison between the constrained experimental and theoretical centerlines was performed using an iterative closest point (ICP) algorithm, as implemented by [39]. The 2-D point clouds used to perform the ICP were the constrained experimental centerline found using the computer vision analysis and the corresponding constrained theoretical centerline found using the LCP analysis outlined earlier in this article. A discretization of N = 200 links was used for the LCP analysis to compare to the experimental centerline. The theoretical centerline was offset by the outer radius of the experimental actuator (1.5 mm) to properly account for actuator thickness.
To facilitate pairwise comparison correspondence, random points were removed from the experimental point cloud until it contained an equal number of points (200) as the theoretical point cloud.
The ICP algorithm then calculated Euclidean distances to the nearest neighbor in the experimental point cloud for every point in the theoretical point cloud. Next, a least-squares best fit transformation was applied to align the experimental point cloud with the modeled point cloud.
The process of calculating Euclidean distances and best-fit transformations was repeated until the change between iterations in average distance between neighbors fell below a threshold of 1 × 10 −7 or a maximum of 2000 iterations was reached.

A. Baseline Computer Vision Measurement Error
The performance of the computer vision program was characterized using ground truth shapes-3D-printed (Formlabs Form  Fig. 19.

1) Observations: Baseline Computer Vision Measurement Error:
The results show that the computer vision algorithm measured actuators with an average baseline measurement error of 1.72% (SD = 1.44%) and 1.88% (SD = 2.60%) for bend angle and bending section length, respectively. The baseline measurement error was relatively constant for both bend angle and length, although performance appeared best when detecting bend angles near the middle of the range. For the proposed application, the baseline measurement error corresponds to fractions of a degree and submillimeter accuracy for bend angle and length, respectively.

B. How Unconstrained Bend Angle Changes With Pressure
The relationship between unconstrained bend angle versus pressure was critical to compare experimental and theoretical centerlines and interaction forces because unconstrained bend angle was used as an input to the model. The average experimental relationship between the unconstrained bend angle and input pressure for each of the three actuators is shown in Fig. 20. A simple linear model was adopted to approximate the empirical relationship between each unconstrained bend angle, θ (degrees), and pressure, P (kPa), for each actuator, described by   Average force behavior displays a knee point at a specific pressure that is different for each actuator, followed by a steeper increase in force as pressure is increased. The knee points were identified using angle-based knee detection. N = 9 samples for each data point. The error bars indicate the standard deviation of the measurements for each data point. process described in Section IV-A. As such, the experimental results were evaluated based on the actuator-specific bend angle to pressure relationship outlined in Fig. 20.

C. Determination of Pressure Range for Comparison
To ensure the middle bending section of the experimental actuator matched the conditions of the model (i.e., symmetry in inflection points created from the first and third bending sections) and that the actuator was in contact with the environmental constraint, only a select range of pressures for each constraint height was analyzed.
The pressure range was determined by plotting the average force versus pressure for each constraint height and identifying the knee point in the force behavior using [40] and [41] (see  120.66-172.37 kPa (17.5-25 psi) for constraint heights of 4.5, 5.5, and 6.5 mm, respectively.

1) Observations: Determination of Pressure Range for Comparison:
The appearance of knee points in the data was most likely attributed to all three bending segments of the actuator not being fully in contact with the environmental constraints. This would explain why the knee points appear at increasing values of the pressure as the constraint height is increased.
As we saw in the theoretical shape analysis as a function of bend angle, we did not anticipate any force being detected for low pressures (i.e., below the knee point). The nonzero forces recorded for pressures below the knee point were likely caused by manufacturing imperfections or gravity causing the actuator to contact the load cell. Although gravity was modeled in our formulation of the problem, we only modeled the middle bending section of the overall experimental actuator comprised of three sections. Because of this, the start and end points of the middle segment may have deviated from the conditions of the model. As such, we only evaluated our experimental results for pressures above the knee point for each constraint height to compare conditions that most closely aligned with the conditions presented in the LCP model.

D. Comparison of Constrained Centerline Shapes
ICP analysis was performed on each image for the pressure ranges determined in Section V-C. An illustrative sample for the constrained actuators is shown in Fig. 22. The illustrative sample was chosen to be the median of all data points (all constraint heights and corresponding pressure ranges) for Actuator 1, with a 5.5-mm constraint height and pressure of 137.9 kPa (20 psi). The results of the ICP analysis are shown in Fig. 23 for both the illustrative sample and all 144 samples.

1) Observations: Comparison of Constrained Centerline Shapes:
The results show that after ICP, the comparison between the theoretical and experimental centerline produced a mean point pair distance and point pair distance error for all 144 samples of 0.31 mm (SD = 0.22 mm) and 1.06% (SD = 0.74%), respectively. The outliers appeared to be caused by the difference in length between the model and experimental centerlines, where the model centerline was consistently longer than the detected experimental centerline. This can be seen in Fig. 23(b), where the point-pair percent error was largest for point pair indices corresponding to the ends of the actuator.

E. Comparison of Interaction Forces
The interaction force for each experimental data sample was compared to the theoretical interaction force from the model. The model used the average experimental bending stiffness of EI = 32.01 Nmm 2 and weight per unit length of w = 2.45 × 10 −4 N/mm found from the free-fold tests. To properly show the force range created by variations in material properties during the manufacturing process, the model was also analyzed using both the overall minimum and maximum bending stiffness and weight per unit length values found from the fold tests. The minimum values were of EI = 15.18 Nmm 2 and w = 2.37 × 10 −4 N/mm, respectively, and the maximum values were EI = 41.80 Nmm 2 and w = 2.47 × 10 −4 N/mm, respectively.
The model also used the average bending section length, , and the average relationship between unconstrained bend angle, θ, and input pressure found from the computer vision algorithm to ensure an accurate comparison between experiment and theoretical models. The results are shown in Fig. 24.

1) Observations: Comparison of Interaction Forces:
The results show that the error in force between the experimental and theoretical actuators trended toward zero as pressure was increased, where the average magnitude of error was 45.34% and 20.98% for pressures at the low and high ends of the tested range, respectively. One possible reason for this is that the interaction forces between the actuator and its environment were minimal at low pressures (i.e., on the order of 1 gf), as shown in the top plots of Fig. 24, so the percent error was large at low pressures. Another potential cause of the trend is that the experimental actuators were not perfectly manufactured, nor perfectly symmetric, so each of the three bending sections in the experimental actuator were not all in contact for low pressures. This asymmetry may have cause the experimental "support conditions" to stray from the simply supported conditions of the model and prevented the end points from lying perfectly along the midline between the environmental constraints. As pressure increased and each bending section made contact with the environmental constraint, the experimental actuator became more symmetric and better aligned with the model assumptions for the support conditions.
These results also highlight the nontrivial effect that small variations in material properties, caused by lot to lot variability in the fabrication method, have on actuator behavior. Although manufactured using the same process and materials, the three actuators varied slightly from each other in structure. This caused differing relationships between both bend angle and force as functions of pressure.

VI. CONCLUSION
In this article, we presented a model for soft actuator behavior in constrained, rigid environments that requires no embedded sensor integration and can be done using only static images with computation times on the order of tens of milliseconds.
We developed a model for soft actuators that deform due to eigenstrains and showed that linear complementarity and Hencky bar-chain methods can be applied to soft robots to provide an elegant and fast solution to modeling in soft actuators. This provided the soft robotics community with a contact modeling method without the need for computationally expensive analysis, such as finite-element analysis (FEA).
We presented a custom experimental setup and computer vision program used to validate our model for both shape and force. This contribution allowed for efficient experimental analysis of a large number of static images. The setup was used to analyze over 297 images (54 ground truth, 99 unconstrained, and 144 constrained images). We showed we were able to analyze images with error below 2% for the critical measurement characteristics. Our analysis produced a mean error in shape between the model and experiment of 1.06% using ICP analysis. Our results also indicated the error in interaction forces between the model and experiment decreased as pressure increased, where the average magnitude of the error was 45.34% and 20.98% for pressures at the low and high ends of the tested range, respectively. This was possibly because the measured forces were less than 20 gf.
Overall, our approach provided speed (∼10-ms compute time) via lower computational complexity at a significant cost of accuracy between displacement and interaction forces outputs. Errors in shape and displacement were dwarfed by errors in interaction force by an order of magnitude. This may be unsurprising as our model and image-based centerline shape acquisition emphasized displacements and not forces. Force measurements were neither made nor required in this method, but predicted shapes match interaction constraint experiments remarkably well. FEA can avoid such an imbalance between force and displacement accuracy, though typically at the cost of more computational and modeling complexity.
Although the experiments of this work focus on hydraulically actuated soft robots, the impact of this work is not limited to the fluid-powered subset of soft actuators. The presented model can be used for any high aspect ratio ( φ) soft actuator that undergoes deformation caused by eigenstrains, independent of actuation technique (e.g., pneumatic, thermal, electrohydraulic, etc.). Overall, this work can be used to efficiently search the design space of soft actuators to determine how actuator shapes and forces are affected by design changes (e.g., material properties, geometries, etc.).
Future work includes expanding this model to account for elastic environments like tissue, as well as an incorporation of a friction model to better capture the environments in which soft robots are often used. Additionally, although FEA is out of scope for this work, future work will focus on a comparison between FEA and the proposed model to further compare the tradeoffs of both methods.

1) Single-Point Approximation With Nonlinear System of Equations:
The gradient descent optimization of a nonlinear system of equations approach is closely aligned with the LCP formulation developed in Section III-C, but has one fundamental difference. Instead of solving for a vector of reaction forces as in the LCP formulation, we instead solve only for the reaction force at the actuator's midspan to minimize computational complexity. This is analogous to using a single "haptic interaction point" to model more complex haptic interactions, as shown by Wei et al. [42]. We specify that the deflection of the mid-span node in the y-direction must be equal to the environmental constraint (see Fig. 25). If the actuator is discretized into an odd number of links, reaction forces develop at the two nodes on either side of the actuator's midspan and their sum is the total reaction force.
We enforce this requirement with the addition of two constraining equations. The first of which is an equation specifying that the sum of all joint angles, φ 1:N , must be zero. This forces both end nodes to lie along y = 0.
The second constraining equation specifies that the deflection of the midspan node must be equal to the constraint height.
In total, the solution requires solving 3N + 2 equations simultaneously (see Table I). The computation times required for both the nonlinear gradient descent and LCP approaches are compared in Table II. A comparison of the percent difference in total reaction force versus bend angle for the LCP and nonlinear system of equations approach is shown in Fig. 26.
A comparison of the actuator shapes generated for the two approaches is shown in Fig. 27.
2) Observations: Comparison of Two Contact Model Approaches: With regard to computation time, we observed from   Tables I and II that the large number of equations being optimized in the nonlinear gradient descent approach caused an increase in computation time that is of order O(10 3 ) to O(10 4 ) larger than the LCP CPU time. The benefit of the LCP approach was that we were able to solve a nonlinear problem using an iterative, linear formulation, which played a significant role in decreasing the computation time of the contact model. The results also showed that the shapes and reaction forces found using the two models are closely aligned, with the disagreement in total reaction force consistently below 0.1% for constraint heights of 4 mm and above. The largest percent differences in reaction force appeared when the constraint height was 3 mm and below. This difference was likely attributed to the noticeable difference in the constrained actuator shape for small constraint heights, as shown in Fig. 27.
Similar to the force comparison, the two methods were closely aligned with regard to actuator shape for constraint heights above 3 mm [see Fig. 27(a)]. However, the performance of the nonlinear gradient descent approach broke down when the constraint height was below 3 mm. We can see in Fig. 27(b) that the actuator shape was not fully constrained within the environment and bowed out on either side of the midspan node. Conversely ,  TABLE II  COMPARISON OF THE NORMALIZED, ACTUAL COMPUTATION TIME FOR THE NONLINEAR SYSTEM  the LCP approach was fully constrained and develops reaction forces in regions on either end of the contact area.