A Lumped-Mass Model for Large Deformation Continuum Surfaces Actuated by Continuum Robotic Arms

Currently, ﬂ exible surfaces enabled to be actuated by robotic arms are experiencing high interest and demand for robotic applications in various areas such as healthcare, automotive, aerospace, and manufacturing. However, their design and control thus far has largely been based on “ trial and error ” methods requiring multiple trials and/or high levels of user specialization. Robust methods to realize ﬂ exible surfaces with the ability to deform into large curvatures therefore require a reliable, validated model that takes into account many physical and mechanical properties including elasticity, material characteristics, gravity, external forces, and thickness shear effects. The derivation of such a model would then enable the further development of predictive-based control methods for ﬂ exible robotic surfaces. This paper presents a lumped-mass model for ﬂ exible surfaces undergoing large deformation due to actuation by continuum robotic arms. The resulting model includes mechanical and physical properties for both the surface and actuation elements to predict deformation in multiple curvature directions and actuation con ﬁ gurations. The model is validated against an experimental system where measured displacements between the experimental and modeling results showed considerable agreement with a mean error magnitude of about 1% of the length of the surface at the ﬁ nal deformed shapes. [DOI: 10.1115/1.4045037]


Introduction
Design and development of smooth, continuous-bodied (continuum) robotics are increasingly aimed at various engineering fields ranging from bio-inspired robotics to medical and healthcare procedures [1][2][3]. Continuum robots, particularly inspired by biology, have become an active research area. Numerous continuum robot arms directly actuated by pneumatic artificial muscles or remotely actuated via concentric tube or and/or tendon-based structures are well established [1,4]. Further studies have focused on kinematics, dynamics, or control of continuum arms and manipulators with results indicating that complex motions are achievable, which can be utilized in a wide range of industrial and healthcare tasks [5][6][7][8][9]. However, the continuum robotics field is no longer restricted merely to the actuation of arms developing a curve or line in space. Rather the application of continuum robotic elements can be extended to actuate spatial surfaces featuring a high degree of flexibility, i.e., so-called large deformation continuum surfaces (LDCSs).
LDCSs have the potential to be widely utilized across a range of engineering applications such as manufacturing, e.g., for providing reconfigurable molds which are currently subjected to costly processes [10,11]. A practical application of actuated surfaces as reconfigurable molds has been previously reported by Habibi et al. [10] in which vacuum-jammed surfaces integrated with pneumatic artificial muscles are characterized to enable reshaping different molds of complex geometries though lacking a reliable kinematic model. Experimental applications such as this require high-resolution models to form the basis of model-based control of the surfaces.
Another application that is currently in use is healthcare tooling to assist patients with the lack of mobility such as soft/flexible exoskeleton systems [12]. Such surfaces also have the potential to be used in aerospace and automotive industries to adapt to and control aerodynamic forces. For all of these applications, the actuation and deformation of the surface must be highly predictable to achieve on-demand, desired profiles consisting of simple to multiple curvatures as schematically shown in Fig. 1.
However, LDCSs have been so far operated mainly based on user intuition and personal expertise rather than on model-based control and simulation. This type of operation would then lead to trial-and-error based methods in their design which raises production costs and degrades performance. To obviate this issue, a flexible and computationally efficient model needs to be developed to accurately characterize configurations due to interactive forces applied by actuators and external force elements. This will enable surfaces to be accurately simulated and the resulting models further made available for model-based control methods. One challenge to model these surfaces analytically is that most of the available beam, shell and plate theories are only applicable when the body thickness (relative to the planar dimensions), and consequently shear effects, are assumed small or zero. Well-known examples include Euler-Bernoulli beam theory and Kirchhoff-Love plate theory which have been developed for thin beams and plates [13]. Although other analytic solutions have been recently developed, so far they are only applicable to 1D structures such as cantilever soft arms undergoing large shear deformation that are subjected to external loadings (e.g., Ref. [14]).
On the other hand, some developed theories for thick-walled beams and plates that account for shear deformations and rotational bending effects, such as Timoshenko beam theory [15,16] and the Mindlin-Reissner theory of plates [17], cannot be used for the large deflection and flexibility considered in this work which are caused by embedded actuators due to linear elastic constraints on their strain-displacement relations [13]. As a result, LDCSs have been solved for numerically rather than analytically. In summary, no basic shell or plate theory has been reported in the literature to be appropriate and applicable for modeling LDCSs that are actuated by continuum robotic arms.
Models proposed for the actuated LCDSs undergoing large deformations include work by Kano et al. [18] which presents a model for a two-dimensional sheet-like robot inspired by the control scheme of the scaffold-based locomotion of snakes. The surface can be bent into relatively large curvatures. However, gravity is not taken into account in this model. In addition, the sheet's thickness and the moment of inertia effects along two coordinates are considered zero. In another work carried out by Medina et al. [19], Euler-Lagrange methods are used for modeling a planar 8 × 8 multilink grid with massless segments. The modeled surface is assumed very thin and highly constrained. The surface model could deform to a curved shape by embedded actuator segments but it cannot tolerate the application of significant forces. The simulation average accuracy reported in the work is 0.25 of a link length.
Another approach called the phantom muscle method was developed by Merino et al. [20] to present a kinematic model for LDCSs when deformed by actuators attached to their edges. This technique uses an infinite number of interpolated curves parallel to the attached actuators. Although this approach introduces a relatively simple mathematical model with good computational efficiency, it lacks the inclusion of several crucial surface parameters such as material properties and gravitational effects.
The work presented here details the development of a kinematic model for LDCSs using a lumped-mass technique that encompasses essential factors such as material properties of a soft surface, inertia forces, gravitational effects, material damping, and in-depth shear effects across thick plates. The use of lumped-mass models has been well reported and acknowledged for their adaptability with large deflections and flexibility [21][22][23], ease of implementation [24,25], and considerable computational efficiency as well as the reliable capture of systems' dynamics [22][23][24][25][26].
In principle, the lumped-mass approach is a close relative to finite element Analysis (FEA) method. Both approaches discretize a body into a finite number of elements which are connected through nodes. However, they are slightly different in the way they treat the nodes and calculate the displacements, which has been well discussed in the literature [26,27].
The model developed here introduces a new application of the lumped-mass approach which will be detailed in Sec 2. Although the approach has shown conformity with large deflections and flexibility well in previous studies [23][24][25], it has never been used as a 3D, two-layer plate model to simulate thick surfaces that include thickness shear effects. Moreover, the integration of such a 3D, lumped-mass LDCS model with a continuum arm model to be actuated in large deflections has not been studied or reported in the literature. The model will first be compared with analytical solutions resulting from classical beam and plate theories when the LDCS is statically deflected under its own weight. An actuator model developed using the Euler-Lagrange method is then integrated with the surface model so that on-demand bending, physical characteristics and geometry of the surface can be evaluated after actuated deformations. It will be shown that the surface model is capable of being bent smoothly into desired profiles for multiple actuator configurations including a single actuator linked to one edge and a pair of actuators linked to two parallel edges of the surface.
In the following (Sec. 2), details of the approach including the model configuration and characteristics, its corresponding equations and a simple theoretical verification of the lumped-mass model will be presented. Section 3 then demonstrates how to model an adaptable continuum robotic arm and link it to the LDCS. Simulated results are then presented for desired bending and actuation profiles of the LDCS. Section 4 provides details on an experimental actuated LDCS setup, with experimental validation of the model for actuation configurations of a single actuator linked to one edge and a pair of actuators linked to two parallel edges of the surface. This section also presents further validation of the model by applying an external concentrated force on the surface while monitoring its transient dynamic performance. Finally, the study will be summarized by concluding remarks in Sec. 5.

Modeling Approach
In Sec. 2.1, the 3D model developed in the current work is explained in detail.

Model Configuration and Characteristics
. The flexible surface model is composed of two interconnected layers; each including a lattice of lumped masses linked together through linear springs. The model also includes springs in all locally diagonal directions around every mass in the model to take shear effects into account and hence represents a more realistic performance of a thick plate or shell than other surface lumped-mass models presented to date, e.g., in Refs. [21,25], wherein only one layer of simple lumped-mass grids has been created. Figure 2(a) shows a representation of the developed model with a close-up of a corner of the surface to display all the spring connections between eight typical masses in the layers. As a result of both the anticipated large deformations and plate/shell thickness, the planes normal to the neutral plane of the surface will not remain perpendicular after the surface undergoes deflection. The resulting model configuration, therefore, links a typical mass in the middle of any layer to eight surrounding masses in the same layer along with nine other masses on the opposite layer of the model as shown in Fig. 2(b). It is not shown in the diagrams for clearness and simplicity, but every line connecting the masses (shown as spheres) is composed of both spring and damping elements (Fig. 2(c)) to absorb energy and provide stiffness elements. The dynamic behavior of each mass in the model is then the result of forces applied by surrounding springs as well as the gravitational force acting permanently in the Z direction. The model can also include the effects of external forces on the individual masses and those at the boundaries with applied constraints, e.g., those with support reaction loads, utilize this external force for their boundary conditions. Likewise, external interactive loads, such as those imposed by connecting actuation elements and/or external environmental forces, are also applied.
There are a few simplifying assumptions in this work. The surface is assumed flat before deformation and its main axes are aligned straight before loading. The rectangular cross-sectional area along any Cartesian direction in the surface is constant before deformation. These sections are originally perpendicular to the longitudinal axes of the surface but they may not remain normal due to inclusion of shear effects and hereupon allowing for motions of different points in different directions. Hence, the in-depth thickness does not need to be considered small. The mass of the surface is evenly distributed between the two layers on the top and bottom of the neutral plane, and the material of the model is considered elastic, homogenous, and isotropic. The dominant deflection of the surface in this particular work is bending about the y-axis in the XY plane, but other motions are achievable. The internal strain energy of a surface segment could be due to bending moment, transverse or axial deflections, and the model is capable of accounting for these distortions.
The equations of motion to determine the dynamic behavior of the flexible LDCS are derived through general Newtonian principles as follows: where U represents the vector of displacement, [M ] the diagonal lumped-mass matrix, [C ] the damping matrix, [K] the stiffness matrix, and F ext and W are the vectors representing external forces and gravity (weight), respectively, acting on the model. To implement Eq. (1), the displacement of masses should be first expressed in terms of the model physical parameters. For this reason, a computation scheme was devised to apply Newton's second law of motion and Hooke's law directly on each mass projected in all the 3D Cartesian coordinates. Figure 3 shows a representative pair of masses, m 1 and m 2 , connected by a linear spring of unloaded length L that move from their origins at time t 1 to other points in the space at time t 2 through an arbitrary deformation of the surface.
The transition of the two masses in Fig. 3 develops a spring force that is actually a component of the term [K]{U} in Eq. (1). This force, according to Hooke's law, is determined as follows: where f q is the force between m 1 and m 2 in Fig. 3 applied through compressing or stretching the connecting spring. Also, q represents q =ê x ,ê y ,ê z that are unit vectors for the Cartesian coordinates X, Y, and Z, respectively. The unloaded length of springs in all the three coordinates q is equally set as l. Likewise, k and c are the stiffness of the spring and damping coefficient, respectively, and u 12 q represents the net stretch or compression of the spring along the coordinate q emerging after deformation, which is determined through the following equation: where u 1 and u 2 are the displacement vectors of m 1 and m 2 defined as u 1 = u 1xêx + u 1yêy + u 1zêz and u 2 = u 2xêx + u 2yêy + u 2zêz , respectively. Hence, u iq (i = 1, 2) in Eq. (3) are the components of vectors u 1 and u 2 in the direction q. Moreover,u 12 q in Eq. (2) is the time derivative of Eq. (3) defined asu 12 q = d dt (u 12 q ). Then, a set of equations whose number depends on the number of lumped masses multiplied by 3 (the number of coordinates) can be solved numerically which was implemented here through the software MATLAB R2016A.
As depicted in Fig. 2(b), a typical mass indexed by "i,j" located in the center of the bottom layer of the surface is bound to other masses by 17 linear springs. The surrounding masses are also denoted by "i, j" formats. The amount of displacement for the central point (i,j) is then given by where δU q is determined through the following term; where the terms u q for the masses in all positions ("i − 1, j", "i + 1, j",…) except the central mass (i,j) are determined using Eq. (3), here noting that the lengths l and L can be assigned different values depending on the position of each mass relative to the mass i,j. In Eq. (5), the index T indicates the masses located in the top layer. Any external load in the model is represented by F (ext)q but the weight of each point mass is denoted by W (= mg) as it is applied only in the direction Z.
It should be noted that the parameters k and c in Eq. (4) have been presumed constant and equal for all springs and dampers between any two masses in any direction, whether axial or diagonal, throughout the LDCS model to make the surface stiffness uniform. In other words, the surface is assumed to be isotropic and homogeneous, which is also the case for the experimental surface detailed in Sec. 4.1. Note also that the terms corresponding to damping forces vanish under steady-state conditions and they affect the displacements only in transition states with no influence on the steady-state accuracy of the results.
It is also notable that the motion described by Eq. (4) should be applied for all the surrounding masses shown in Fig. 2(b) and consequently for all the existing masses in the lumped model to be solved at each time iteration. The equations for the masses located at the top layer are similar to Eq. (4) where the places of the top and bottom layers are reversed. Here, the index T is replaced by B indicating the masses located in the bottom layer of the model. Also, the signs +/− in some of the existing terms are changed accordingly based on predefined geometrical positions. It is worth mentioning that the equations of motion for boundary masses vary based on imposed conditions. When a mass is constrained in one direction, the corresponding displacements must become zero in that direction rather than being calculated through Eq. (3). Likewise, for the masses directly linked to actuators, the amount of displacement is initially dictated by the motion of an actuator to which the mass is bound. To clarify this issue, the development of an actuator model by which the LDCS is bent into the shapes desired for the present work is presented in Sec. 3.1.

Theoretical Model Verification.
Having developed the surface model, the LDCS can now be actuated through specific configurations using a continuum robot arm attached to the surface. However, prior to this, the surface-based model was tested, loaded under its own weight, to evaluate its consistency with results yielded by Timoshenko beam theory. The 3D lumped model shown in Fig. 4(a) is composed of 24 lumped masses, total length L = 0.5 m, width w = 0.1 m, cross-sectional A = 0.01 m 2 , total mass m = 0.3 kg, Poison's ratio ν = 0.5, and Young's modulus E = 40 kPa deflecting under its own weight. Given these values, using the approach presented in Ref. [24], developed for planar lumpedmass arrays, along with using equivalent spring constant for series/ parallel springs, an average value of k = 250 N/m was worked out as the spring stiffness matching the properties of this model.
Note that not all the spring-damper links are shown in the figures pertinent to the lumped models in this work. In the presented Fig. 4(a), only the links on surface and edges are displayed for better visibility and to avoid complicating the images.
The results of this initial test, depicted in Fig. 4(b), highlight a reliable conformity between the developed lumped-mass configuration of the surface and the Timoshenko theory, with a maximum error of 1.1 mm at the end point, in comparison to the model dimensions of 500 × 100 × 100 mm. Note that the maximum deflection of the beam calculated and shown here is less than 10% of its total length. However, when the weight of the model was increased to cause larger deflection and curvature, the two result sets found further departure from each other due to the inconsistency of Timoshenko theory with such large deformation as mentioned in Sec. 1. Hence, a different numerical method or an experimental test was needed to validate the model when undergoing larger deformations.

Simulation of Actuated LDCS
This section describes how the modeled LDCS is bent and deformed into the profiles of interest followed by the presentation of initial simulated results after actuation. 3.1 Integrated LDCS-Arm Model. Deformation of the surface in this work is provided by a controllable continuum arm linked to the surface and used as an actuator. Here, it is not our aim to present a novel development in modeling such actuators, rather we illustrate how the surface model is moved and deflected into the shapes of interest with certain curvatures to be compared with the empirical tests accordingly. The development of soft actuators themselves is a crucial phase in designing robotic LDCS surfaces [28]. An actuator model was developed here to match the physical realization of the surface model. The general concept of the actuator model is shown in Fig. 5(a), and this is then embedded into the surface model. The general idea of this arm model has been previously developed and validated by other researchers, e.g., a study carried out by Giri and Walker [29] wherein a section of a continuum arm is modeled using lumped model elements and application of Lagrangian principles. The main difference with the presented model is that the arm sections in Ref. [29] are driven by input forces representing air muscles in the system while the current work applies torque joints between the lumped segments to rotate them and bend the entire arm resulting in fewer actuated degrees of freedom.
The central backbone of the model is assumed able to rotate along the segments for generating smooth bends. As a result, the two marginal, parallel edges of the model surface can elongate and/or contract while bending. These two edges will be then linked to the top and bottom borderlines of the surface model.
The actuator model is a multi-body system consisting of lumped segments of mass m l and moment of inertia I as shown in Fig. 5(a). The segments of the core spine are joined together via torsional springs and torsional dampers whose stiffness and damping coefficients are denoted by k t and c t , respectively. A clearer configuration of the model that incorporates all necessary characteristics for formulation and matches well with the modeled LDCS is the 4-link, lumped parameter model illustrated in Fig. 5(b) where its flexibility is due to the fitted torsional springs. In our case, the arm is restricted to planar motions, similar to the most of previously developed soft bending actuators (e.g., Ref. [30]), which results in movement in the XZ plane, as this is also how the experimental arm is setup. The transverse links A 1 B 1 , A 2 B 2 , A 3 B 3 , and A 4 B 4 are attached and perpendicular to the links L 1 , L 2 , L 3 , and L 4 , respectively, and hence their motions determined by the corresponding links. In this model, the length of all the related links is assumed equal as L i = l l and |A i B i | = l b .
The major variable of this system is θ i , the rotation angle of each link in the XZ plane, through which the state of the system can be fully described given that the left-hand node of the first main segment (O 0 ) and the transverse link A 0 B 0 are fixed in space. As both the kinetic energy (E k ) and the potential energy of the actuator (E p ) can be evaluated with respect to θ i , its motion can be described through a Lagrangian formulation via the following: where L n (= E k -E p ) is the composite energy term and T r denotes the external torque acting on the coordinate θ r . As shown in Fig. 5(b), the torques (T r = T 1 , T 2 , T 3 , and T 4 ) are applied on the lefthand joints of the links at O 0 , O 1 , O 2 , and O 3 , respectively, to generate a desired curvature throughout the arm. Having calculated E k and E p , the Lagrangian term for the four-link arm shown in Fig. 5(b) is given by + 5θ 1θ2 cos(θ 1 − θ 2 ) To evaluate the equations of motion for this system, the masses of the links were considered in the potential energy and the damping effects were applied through the principle of virtual work and included in T r . Given the energy term in Eq. (7) and substituting in Eq. (6), it is then possible to derive the equations of motion for the arm which are solved numerically using MATLAB. Then, the rotations of all the four main links (θ i ) and, consequently, the displacements of the transverse links i.e., the points O i , A i and B i shown in Fig. 5(b) are determined using geometrical relations. For example, the XZ coordinates of the point A 3 are given by x A3 = l l (cos θ 1 + cos θ 2 + cos θ 3 ) − l b 2 sin θ 3 z A3 = l l (sin θ 1 + sin θ 2 + sin θ 3 ) + l b 2 cos θ 3 (8) which are then subtracted from their original values to arrive at net displacements. Subsequently, the points A i and B i are connected to the corresponding lumped masses located in the top and bottom layers of the surface model, respectively, as illustrated in Fig. 6. This connection is made through stiffer springs (as the values are given in Sec. 3.3) to enable force-based connection between the elements. Note that every three pairs of the lumped masses of the surface in the longitudinal direction (X ) are bound to one main link (segment) of the actuator as seen in Fig. 6. The intervening masses are not joined to the arm. As a result, the surface and the continuum arm are bound together at five sections indicated by Figs. 5 and 6. It is mentionable that the arm model was initially developed with the same number of masses along the surface edge to connect every link's end to a mass on the surface model. However, that arrangement made only a small improvement on the accuracy of results, less than 0.2%, at the cost of a significantly increased runtime. Because the purpose of this model is for use in model-based control, computational efficiency applications, reducing computation time over absolute accuracy is an essential factor. To match the two models, the thickness of the surface should be adjusted equally with the length of transverse links, i.e., |A i B i | = l b = l z . Moreover, since the length of each link of the actuator is twice the distance between every two masses of the surface in the x-direction; in all simulation results, it has been assumed that L i = l l = 2l.

Simulation
Results. The dynamics of the developed flexible grid when actuated is influenced by gravity and the external loading applied by the actuators integrated with the surface. In this work, two configurations for the integrated LDCS-arm model were developed and tested: Test 1-one actuator is mounted on one edge of the surface (Figs. 7(a) and 7(b)) and Test 2-a pair of the actuator model are linked to two parallel edges of the surface (Figs. 8(a) and 8(b)). In both arrangements, the surface is clamped along the indicated edge to hold the surface up from the ground.
The surface model contains a two-layer square grid of 9 × 9 masses developed according to the modeling approach described in Sec. 2. Thus, it consists of 162 masses of m = 1.5 g with the total mass M = 243 g. This conforms to the characteristics of experimental test surfaces manufactured for the test setup. The total length and width of L tot = 0.160 m (i.e., l = 0.020 m), thickness l z = 0.010 m, and, therefore, a cross-sectional area A = 16 * 10 −4 m 2 were selected. A 9 × 9 grid was chosen to represent the surface as an initial trade-off between accuracy and computational efficiency of the system as discussed in Sec. 4.3.
A value of k = 400 N/m was determined as the spring stiffness of the LDCS model to match the Young's modulus and Poisson's ratio of the flexible surface in the experimental test described in Sec. 4.1 using the approach presented in Ref. [24], previously developed for planar lumped-mass arrays. The value for the springs linking the arm(s) to the surface was selected as k a = 4000 N/m to represent a firm connection between the two models.
Physical properties of the continuum arm model were selected to match the characteristics of the fabricated actuator that will be discussed further in Sec. 4.1. Consequently, each link of the model was given a mass of m l = 0.01 kg and length l l = 2l = 0.040 m, while the massless transverse links (A i B i ) were fitted equal to the surface's thickness, i.e., l b = l z = 0.01 m.  As can be seen, the top center of the surface has slightly sagged down under its weight which is due to symmetry in the boundary conditions applied by the two actuating arms. In both Figs. 7 and 8, the actuating arm is shown in black lines in the lateral side(s) of the surface model where it has displaced the surface from its original position (as indicated in the figures) placed on the horizontal XY plane and bend it up to obtain the curvature determined by the values given to the applied torques (T i ) and the stiffness of torsional springs (k ti ). Note that since the simulation is initialized from zero gravity, the actuator model does not need to be pre-strained and pressurized to keep the surface straight in the XY plane. In other words, actuating and gravity forces are applied simultaneously at the beginning. It should also be noted that Figs. 7 and 8 have been drawn on purpose slightly detached to highlight the distinction between the two models, when in the model they are actually connected by the spring k a .
Further details on the testing procedure are given in Sec. 4 where the equivalent experimental tests are presented to provide ground truth for evaluating the modeling results.

Test rig.
To validate the modeling results, a test rig, shown in Fig. 9, was set up consisting of an aluminum frame, position sensor system, a personal computer data acquisition (PC DAQ), and a pneumatic system to operate the actuators. The sensor system is a 3D Guidance trakSTAR (Ascension Corp., Shelburne City, VT) chosen and used for its convenience and high accuracy with measuring displacements up to 0.1 mm. The sensor system's main box was connected to the PC using a USB cable and three sensor probes (Model 180 were plugged to the main box. The sensors can be inserted in the surface at any desired position. The sensor system also includes an electromagnetic generator (MRT) to generate a magnetic field working in conjunction with a transmitter that establishes the coordinate frame and tracking volume. The system is then capable to sense the displacement of every probe in the magnetic field in three spatial coordinates. A soft surface was then fabricated from Ecoflex-0050 silicon rubber poured into a rectangular mold to be solidified and shaped into the same effective dimensions considered for the developed LDCS model. Also, a clamping holder to fix one edge of the surface was built through prototype 3D printing. Similarly, the continuum robotic arm actuator was made from a material known as DragonSkin-0030 silicon rubber formed through a designed mold shown in Fig. 10(a). Manufacturing the arm was then accomplished by fitting reinforcement fibers and an inextensible layer to be adapted and integrated with the fabricated surface as presented in Fig. 10(b). The details on the method to fabricate this type of actuator are provided by Polygerinos et al. [30]. The actuator is hollow and operated by air pressure which enables to bend up in different desired curvatures. The bending curvature is proportional to air pressure and controlled by varying the input analog voltage of a proportional valve (SMC ITV2000) in the pneumatic system. In this experiment, air pressures of 81 kPa and 116 kPa were used, respectively, in Test 1 (to bend the actuator into a curvature of radius r = 0.131 m) and Test 2 (to bend the actuators into a curvature of radius r = 0.092 m). The curvatures were empirically obtained via fitting the position points measured along the length of several points along the curve. This consistency in curvature is also in agreement with the beam theory and modeling method of the soft actuator from Ref. [30]. Once the soft surface is attached to the surface, it can be approximated as a uniform payload acting on the arm causing a constant curvature shown in Figs. 13 and 14.
The tests were carried out so that the manufactured surface and the LDCS simulation undergo the same loading and boundary conditions as well as geometrical and material properties. The primary difference is likely in the material damping of the modeled LDCS, which was considered in this study as c = 1 N s/m to improve performance of the model, effectively by not letting the model vibrate forever, while maintaining results similar to the experimental values. This difference is clear in Fig. 19  Testing and Tension/Force Gauge, China) shown in Fig. 11(a) was utilized. This machine applied a uniform axial force on the surface, gripped at one edge and gently stretched, with applied force and resulting displacement measured. The resulting data were used to draw the force-displacement graph shown in Fig. 11(b) to calculate Young's modulus as E = 75.35 kPa for this case. In addition, the Poison's ratio of the surface, made from the incompressible material, silicon rubber, is considered to approach 0.5 as was also found in Ref. [31].
To do this tensile test, a surface with the same dimensions of the model (160 mm * 160 mm * 10 mm), was used to provide the Young's modulus of the surface. As shown in Fig. 11(b), the strain obtained in this test is then less than 7%. Note that within this range (generally the strain within 0-1), the stress and strain relationship in hyperelastic materials such as this type of silicon rubber is almost linear, and it can be expressed by Young's modulus [32][33][34]. The large deformation experienced by our surface is in the overall shape and spatial movement of the points across the surface that leads to considerable overall relative bending. However, the distance between every two marked points on the surface do not locally experience a large elongation or contraction that would lead to typical large elongation in hyperelastic materials, e.g., a common range between 100% and 400% [35] [1]. Because the elongation of the surface measured in this experiment does not go beyond 10% shown in Figs. 13 and 14, we may consider our surface to have very small elastic strain for its given material [32][33][34].

Testing Through Experimental
Procedure. The physical parameters such as mass, length, and other dimensions in both experimental tests were chosen to be the same as those in the developed models for the simulation test. However, the mechanical properties, i.e., Young's modulus and Poison's ratio, were first found empirically as explained in Sec. 4.1, and then applied to To proceed with the experimental test, once the arms are operated by air pressure, the actuator(s) of the LDCS reshape it into desired curvatures. Since the actuator(s) in the experimental setup operate(s) almost linearly as a function of input pressure against curvature [30], the supplied pressure was increased gradually to reshape the surface into the curvature resulting in the lumped model.
To determine bending level and curvature for the arm models, the torsional springs pinned at the joints shown in Fig. 5 and O 3 ) are given specific values to result in a bend on the arm accordingly due to the equal, constant torques applied on the joints (T 1 , T 2 , T 3 , and T 4 ). In other words, since same amount of torque is applied on these four joints connecting the main links of the arm(s), when a proportional set of spring stiffness are chosen and allocated to the linear torsional springs, each link rotates relative to its previously positioned, adjacent link with an angle exactly proportional to their stiffness difference. For example, this set in the case of Test 1 was selected as k t1 = 12, k t2 = 6, k t3 = 3, and k t4 = 1.5 (N/m), while the applied torques were selected as T 1 = T 2 = T 3 = T 4 = 0.5 (N m) which provided the desired curvature. In fact, this caused each link to rotate as twice as much as the previous link to finally achieve the overall curvature of radius r = 0.131 m (k r = 7.63 m −1 ) for the continuum arm.
To detail the displacement of the LDCS during actuation, six points across the surface were indicated and labeled on its unloaded state as shown in Fig. 12. These points are assumed to be located in the central, neutral plane of the surface which is averaged between the masses positioned in the top and bottom layers. The probes of the sensors used for the measurements are very sharp, slender, flexible and can be inserted into the soft surface directly. Hence, it is assumed that they have very little effect on the mechanical properties of the soft surface.
Due to a limited number of available sensors in the experimental tests, the three sensors in both setups (Test 1 and Test 2) were first attached to the points P 2 , P 4 , and P 6 to measure their displacements. Then, the same sensors were attached to the other three target points, i.e., P 1 , P 3 , and P 5 , to obtain the desired displacement points used for comparison.

Empirical Results and Comparisons.
The results of Test 1 and Test 2 are presented here to evaluate the validity of the model. Figure 13 which is related to Test 1 for the experimental surface, corresponds to Figs. 7(a) and 7(b) of the simulated results. Figure 14 indicates Test 2, corresponding to the simulated results in Figs. 8(a) and 8(b). In both cases, the final surface displacement is shown. Table 1 shows the resulting displacements for the two sets at each point for Test 1.
As points P 1 and P 2 in Test 1 are firmly joined to the continuum arm, and this does not twist in the y-direction, their displacements in the y-direction (u y ) remain zero. However, these points in the simulation do have minor transitions in the y-direction due to the simulated spring connections between the surface and actuator. Figure 15 depicts the data in Table 1 via a bar chart to visualize how closely modeling and experimental results confirm each other in different coordinates.
Similar to Test 1, the displacement data for all six points in Test 2 are shown in Table 2. The results in Test 2 from the two presented sets show that the average absolute error is less than 1 mm. As    before, some points in the y-direction (u y ) are not displaced. This is because the four points P 1 , P 2 , P 5 , and P 6 are attached firmly to the two parallel actuators. Also, the two points P 3 and P 4 are located in the middle of the surface and which are then subjected to the symmetrical boundary and physical conditions about the surface's central X-axis. Due to this symmetry, as expected, some identical displacements for a few pairs of points result in Test 2 as can be seen in Table 2 Similarly, the data in Table 2 have been utilized in Fig. 16 to plot a clearer picture of agreement between modeling and experimental results obtained in Test 2. Table 3 shows a summary of absolute error between modeling and experimental results at six points of the surface considering one decimal digit for both of the implemented tests. Shown in Table 3, apart from the absolute errors in each coordinate (u x , u y , and u z ), are the spatial magnitude of these displacement errors for each point (δ u ) determined through δ u = u 2 x + u 2 y + u 2 z . As seen in Table 3, the measured points in Test 1 results in very large deformations and therefore reveals accordingly larger errors (mean absolute error of 1.3 mm among the coordinates u x , u y , and u z ) than Test 2 (mean absolute error of 0.7 mm among the coordinates u x , u y , and u z ). Overall, the average absolute error for both tests was found to be less than 1 mm among the coordinates u x , u y , and u z , which is very small in comparison to the overall surface dimensions of 160 × 160 × 10 mm (less than 1% of the length of surface's sides). In addition, the maximum spatial magnitude of the errors (δ u ) is 2.8 mm that corresponds to point P 6 in Test 1. However, the mean value of δ u from the both tests was found as δ umean = 1.75 mm which is still considered small as it is only slightly above 1% of the length of surface's side (1.06%).
One main reason for this small error is the limited number of nodes (masses) used in the model here for computational efficiency whereas with other discretization model methods, increasing the number of nodes would lead to more precise results. However, achieving an optimal number of nodes to be used for the LDCS model requires a comprehensive optimization process that is out of the scope in this work.

Further Validation: External
Loading and Dynamic Transient Performance. As mentioned above, one of the main goals for the modeling of the actuated surface is to evaluate against external loading in static and dynamic conditions. For this reason, in addition to the gravitational effects, an additional mass representing a concentrated, constant external force was applied to the top, center of the surface at point P 4 in the configuration of parallel arms as shown in Fig. 17. This image shows the final deformation of the integrated actuator-surface after the transient operation of actuation ended.
The model configuration in this test was again composed of two lattices of 9 × 9 masses or 162 masses in total across the surface. The value of each mass was chosen as m = 9.87 * 10 −4 g, i.e., with the total mass of M = 160 g. The additional mass was selected as m 0 = 40 g to quantify the external load as F 0 = 0.39 N.
The simulation results of this test are presented in Fig. 18 wherein the two-layer LDCS model has been bent up by two parallel continuum arms of identical physical properties embedded in the surface. The actuating arms, shown as black lines in the lateral side(s) of the LDCS model, have actuated it into a curvature of radius r = 0.115 m (k r = 8.69 m −1 ) at the side edges. It can be seen that the top center of the surface where the additional mass is attached sags down due to the concentrated external force applied to this point. In the experiment, to ensure that the actuation and gravity forces operate at the same time, two plates as shown in Fig. 17    the surface at the start. When the actuation force was applied, the surface was then lifted from the support. Similarly, in tests with the additional weight, it was held by a movable seat support to ensure it would not hang off of the surface before actuation, which results in a starting displacement of 0 at t = 0 s, as shown in Fig. 19. The seat support was moved aside at the same time as the actuation started. Figure 19 shows the displacement of point P 4 on the surface versus time in the x-, y-, and z-directions for both the simulated and experimental tests. As shown, the point undergoes tiny fluctuations in displacement (with a maximum vibration amplitude of 3 mm for the 160 mm-long surface) in the initial transient period due to the sudden movement and low material damping (c = 1 N s/m as explained in Sec. 4.1). Differences in this transient region are small and likely due to material property differences for c and k in the model. After the transient dies out, the two results converge and settle very close together in the x-and z-directions, indicating a reliable static performance for the developed model. The results in the y-direction remain zero as expected due to symmetry in geometry, loading, and boundary conditions applied by the two parallel arms positioned equally apart on its two sides.

Conclusions
This paper has introduced and validated a novel 3D, two-layer, lumped mass-spring-damper model to describe the behavior of actuated surfaces undergoing large deformations. The study has also extended the application of a lumped-mass approach for characterizing and representing thick flexible plates in 3D space where a continuum robotic arm and flexible surface are integrated together. The full model takes into account interactive forces (between the actuating arm and surface), as well as physical and mechanical properties of the system such as mass, elasticity characteristics, gravity, material damping, and in-depth shear effects. The static deflection of the developed surface model under its own weight was first compared with the well-known Timoshenko beam theory with a maximum error of 1.1 mm at the end in comparison to the model length of 500 mm. A test rig was constructed for two simulated surface-arm configurations, Test 1-Single actuator along an edge and Test 2-Parallel actuators along two edges, for experimental comparison purposes. In addition, the developed model accounted for deformations resulting from a combination of loads applied by the actuation arms, gravity, and external forces, while still accounting for in-depth bending shear effects of thick flexible plates. The model further successfully demonstrated the transient dynamic performance of actuated surfaces undergoing large deformation while experiencing concentrated external loading.
Simulated and experimental results show that the model is capable of accurately predicting profiles and curvatures of the actuated LDCS due to applied forces by the continuum arm(s), whether dynamically over the transient actuation time or statically after the end of motion with a mean error magnitude of about 1% of full surface length at final deflected positions.
In summary, the proposed model is a new methodology to enable a modeling method for actuated surfaces undergoing large deformations through the use of continuum actuation. This approach is primarily focused to present a middle ground between FEA techniques (usually with an insufficient level of computational efficiency) and very simple analytical models (usually with low level of accuracy or incapable of modeling shear deformations) for such structures for use in model-based control methodologies. In line with this characteristic, future work will focus on the trade-off between computational efficiency and accuracy of modeling when applied to model-based control methodologies, particularly in applications where high computational power is not available. The model presented here could then be used in model-based control strategies across a range of highly deformable continuum robotics applications such as the manipulation of parts in manufacturing environments, soft/flexible exoskeleton systems in health care, and deformable surface control in the aerospace, automotive, energy and food processing industries.   18 Simulation results for the lumped-mass LDCS model undergoing a concentrated external force displayed at its final static deformed shape after actuation by two continuum arms embedded in its two parallel edges Fig. 19 Comparing the dynamic behavior of the modeled LDCS and the empirical surface through recording the displacement of point P 4 as indicated in Fig. 12 overtime in the three directions X, Y, and Z