Symmetric Kinetostatic Behavior From Asymmetric Spatially Curved Beams

Numerous studies have been conducted on the design of complex monolithic spatial elements with the ultimate goal of increasing a structure’s rigidity by reducing elastic deformations [1]. On the other hand, there are various applications where large-deformable elements with specific force-deflection behavior are desirable. Thus, introducing compliance along the part becomes essential. These monolithic compliant elements are classified as a type of compliant mechanism (CM) as they can provide a specific motion while they are elastically deforming [2–4].

Several engineering applications could benefit from monolithic CMs with a spatial range of motion as a more efficient replacement of complex conventional (rigid body chain) mechanisms, applications where slenderness and being lightweight are essential, e.g., gravity balancers for wearable devices like exoskeletons and prostheses, where a certain compensating force within a spatial motion is expected in a very limited design space around the user’s body [5–8].

In recent years, spatial CMs have become a topic of interest in the CM design community as planar CMs have been extensively investigated. Yet more complex design requirements were defined, which could only be addressed by three-dimensional motions. Efforts have been made in both the characterization, as well as in the design and optimization of spatial CMs. For instance, in flexure mechanisms, there are plenty of examples and principles [9]. However, most of the developed spatial mechanisms are for precision devices, which usually consist of several parts that make them bulky compared to their output range of motion.

Designing monolithic spatial CMs with an extensive range of motion was emerging in recent years [10–13]. The main advantage of having a fully compliant monolithic part to obtain a desired kinetostatic behavior is that the whole shape effectively takes part in the CM’s characterization.

There are several widely used CM design methods, e.g., freedom and constraint topology (FACT), building block approach, Pseudo rigid body model (PRBM), etc. [14–16]. However, using these methods for designing large deformable spatial monolithic parts are not feasible, because capturing the distributed contribution to the endpoint stiffness, in addition to a desired kinetostatic behavior from a spatially curved shape is more complex than what these methods can propose. Another approach is to take advantage of automated design methods like optimization. Concerning the complexity of the defined requirements, the latter approach was assumed to be a better choice to find these compliant shapes.

In this regard, there are a handful of studies on the topology optimization of monolithic CMs [17–20]. The output of these works were usually rather complex monolithic shapes which make them only fabricable with additive manufacturing techniques.

To achieve complex kinetostatic behavior and yet keep the CM manufacturable with common techniques, choosing a simple topology and optimizing its shape parameters can be considered a better solution [21]. There is extensive research in this regard, mainly in planar mechanisms [22,23]. However, making use of spatial mechanisms can enable other possibilities, e.g., torsion of elements.

In this regard, spatial geometries varying from rods to shells were scarcely investigated [24–28]. There is no doubt that reaching exceptional kinetostatic behavior might be possible with different forms of compliant elements. However, from a utilization perspective, as well as a fabrication perspective, it appears convenient to consider the simplest topology of a cantilevered spatial beam as the monolithic compliant element.

A useful yet undiscovered kinetostatic behavior from this type of CMs is sought in this work, where the endpoint of a spatially curved beam has symmetric stiffness behavior in its plane of motion, even if the beam is asymmetric in the sense that it is not grounded in the symmetry plane of its stiffness field. This situation emerges in several applications of manipulator-type devices where the endpoint symmetry plane is occupied and no room is available for the grounding point on that plane, these devices can be widely used in surgical robots and aerospace structures. Also, in wearable devices where the desired symmetry plane of the limbs’ motion is occupied by the body itself, or when a symmetrical support force is required but a symmetrical structure is not possible. For instance, consider a unilateral upper-body bending support for wheelchair users, it is not safe to lock the user’s upper-body in a wheelchair with a bilateral, symmetrical structure, and a normal unilateral, asymmetrical structure cannot provide the symmetrical support on the chest and will turn open upon forward bending. In this instance, an asymmetric beam with symmetric kinetostatic behavior, such as the one proposed in this paper, can supply the required bending support force while not restricting the user’s freedom. The majority of the parameters used in this study are derived from the latter application. In our earlier work [29], we showed that there exist non-uniform distributions of cross-sectional properties along an asymmetric beam for which the natural rotating behavior was corrected toward more symmetric endpoint behavior.

This paper aimed to realize the effect of tuning the shape and sectional properties of a spatially curved asymmetric beam to enhance the symmetry of the endpoint’s stiffness in a specific range of motion. By utilizing B-spline, a parametrization scheme was developed to generate beam shapes in the prescribed design space from a small set of parameters. An objective function was defined to evaluate the symmetric kinetostatic behavior of the endpoint in the range of motion. Additionally, an optimizer was used to improve the objective function iteratively by tuning shape parameters. Finally, experiments were conducted to validate the kinetostatic performance of the shapes generated by this process.

The paper is structured as follows. In Sec. 2, the definition of the requirements is specified together with the description of the optimization process for different combinations of parameters. Also, the formulation of the objective function and the optimization settings and the finite element solver are described. Furthermore, the verification procedure using experiments is discussed. In Sec. 3, the resulting beams from all cases and their performance are shown and compared with the benchmark beam. A discussion on the validity of the results and possible improvements are given in Sec. 4, and the conclusion is given in Sec. 5.

The process of making and evaluating beams has different aspects. These aspects are discussed in the following subsections. The overall procedure can be summarized as follows: (1) A scheme for parametrizing the shape of beams using B-splines was developed to form the beams’ shape using a limited number of control parameters; (2) the beams were generated by optimization of different sets of shape parameters; (3) the objective function of the optimization was defined to assess the symmetric kinetostatic behavior in the desired region; (4) the displacements of the beams under the objective function’s loadings were numerically computed using a self-developed finite element model; and (5) an experimental setup was designed to verify the beams’ performance.

As discussed in Sec. 1, the topology of the simplest form of a monolithic CM is that of a cantilever beam. This includes a fixed clamp in the base and an asymmetric spatially curved beam which ends at an end-effector point at the tip.

The position of the grounding point was assumed to be at the origin of the Cartesian coordinate system, and the beam endpoint was chosen arbitrarily to be at X_e = 0.20 m, Y_e = 0.15 m, and Z_e = 0.50 m as shown in Fig. 1. The symmetry plane was selected as the plane that coincides with the endpoint, parallel to the YZ plane. The beam shape was formed based on the set of parameters described in the following subsections.

In this work, the parameter optimization to form the beam in each of the iterations was divided into two branches, the first branch contains parameters related to the beam’s shape and the second set contains sectional parameters and orientations along the beam. The shape optimization was implemented in three steps: (1) a fixed-shape beam, (2) an optimized free-form shape beam, and (3) an optimized shape under a curved planar constraint.

The second set concerns the cross section and has two possible conditions: (1) a beam with a circular section, which means that only the circular beam’s radius was optimized, and (2) beam with an I-section where all sectional parameters and sectional orientations of control points were optimized, as shown in Fig. 2. The combination of these two sets made 3x2 beams (see Fig. 5). These six beams are described in the following subsections.

It is important to note that, since we used a beam model in the implemented finite element method (FEM), the sectional dimensions were ultimately reduced to four main parameters, i.e., the area, the second moments of inertia about the two main axes, and the torsional constant (A, I yy, I zz, J). Thus, it might ask why we didn’t optimize parameters instead of sectional dimensions. However, these parameters are not independent, and thus, optimizing them directly as separate variables could lead to unfeasible results. For this reason, it was chosen to optimize the dimensional parameters of a prevalent section. Among those, the I-section was selected since changing its dimensions enables a large variety of combinations of these four principal beam parameters. However, as the optimization is essentially about determining a collection of four principal beam parameters (A, I yy, I zz, J), another non-bisymmetric section with a relatively low torsional constant can do the same task as the I-section.

In the second beam, the dimensions and orientations of an I-section along the beam were subjected to the optimization, and the beam shape remained the same as it was described for Beam 1. The I-section was defined by their web height H, flange width B, flange thickness h, web thickness b, and the section orientation φ (see Fig. 2). The first four parameters were subjected to optimization once for the whole beam. The orientations were optimized at five control points, and all the other sectional orientations along the beam were interpolated based on the values of those five points.

The beam shape’s optimization process was based on finding the best positions of intermediate control points to form the beam shape. The shape was made based on a B-spline between the grounding point and the endpoint described in subsection 2.1, concerning the optimized intermediate control points. The positions of all other nodes along the beam were interpolated from this B-spline, and together they made the shape of the beam. Theoretically, it was possible to choose any number of control points to make this spatially curved line. Here, we used five control points in the design space.

Bounds on the search space were applied to avoid the shape crossing itself or going out of the design space. The section of Beam 3 is again a filled circle, with one optimized radius. It is important to note that in all beams, the optimization process for section and shape parameters was performed simultaneously. Here, in Beam 3 a set of 16 parameters, including three coordinates for five control points and one radius were optimized.

This beam was made up of a combination of both the mentioned procedures for shape optimization and sectional parameter optimization, which made it the most advanced beam among all six. These two sets of parameters were optimized together and enable control of all possible shape parameters of this scheme. A matrix consisting of three coordinates for each of five control points and a set of four variables for sectional dimensions and five parameters for orientations of sections were optimized together. These 24 parameters together formed the beam’s shape.

The procedure for this beam was like Beam 4, only with two instead of three coordinates for control points to meet the constraint as defined in Eq. (4). It is important to note that in this beam, not only the control points, but also all other interpolated points along the beam were defined to obey this equation to keep the beam shape always on this extruded ellipse constraint, as shown in Fig. 3.

A degree-four B-spline was used to interpolate the shape and the sectional properties based on the optimized control points. The B-spline was chosen to have an open uniform knot vector to ensure that the first and last optimized control points coincide with the first and last parameters of the beam itself, which were the grounding and the endpoint. The internal knots were determined based on the de Boor algorithm [30] to equalize the weight of all control points. Using B-spline interpolation ensured a smooth transition between optimized parameters of the beam and avoided any discontinuities that could have happened due to sudden dimensional changes in the finite element model.

All the displacements under different loading conditions were derived from a self-developed finite element solver. The developed solver uses geometrically non-linear co-rotational beam elements. The mentioned beam elements were introduced by Battini [31], based on the Euler–Bernoulli formulation. This solver was selected as it can handle non-linearities due to large elastic deformations. Furthermore, having access to the source code of the solver made it an easier choice to tailor it for the aim of this work. For the optimization process of the parameters, the Multi Start option from the matlab^®optimization toolbox was used. Five random starting points for the fmincon function with the Interior-Point algorithm were set. With an objective function tolerance of 1e − 5, the maximum number of iterations is set at 100. Upper and lower bounds were used to keep the sectional parameters in a way to preserve the I-shape of the resulting sections and keep the beam shape inside the design space. The material constants for all beams were Young’s modulus of E = 200 GPa, and shear modulus of G = 76.9 GPa.

For verification of the method’s efficacy, the resulting beams’ kinetostatic behavior has been checked with experiments. To make the results from this verification comparable with the results from the developed code, the same loading conditions from the objective function were implemented on the beams, and the errors were derived from the same objective function of Eq. (1).

All six resulting beams have been 3D-printed with the multi-jet fusion method and polyamide 12 (PA12) as the print material. The beams have been scaled by 0.7 due to limitations in the printing size. To generate endpoint loads from different directions, the beam grounding was set to rotate about a horizontal axis. The load was exerted by known weights at the endpoint, which was always pointing vertically downward. The experimental setup is shown in Fig. 4(a). The resulting displacements have been extracted using vision-based measurement between the loaded and unloaded tip positions.

To evaluate the performance of the described method, the resulting objective functions of all beams with different optimized parameters are shown in Table 2.

The first set of errors directly resulted from the optimization procedure and self-developed finite element model. The second set was based on displacements obtained from the experiments.

The isometric and top views of the resulting beams are shown in Fig. 5. To represent the kinetostatic behavior, in the field of action, an isoforce mapping at the endpoint of the beam is presented with black lines, which represents the displacement of the endpoint under the constant magnitudes of the force F = 20, 40, 60, 80, 100 N and a full cycle of the angle α. The mirrored colored lines represent the endpoint path under increasing force from zero to 100 with mirrored angles ±α shown in the same color. The isoforce mapping together with the colored mirrored loading mapping can indicate the symmetry in the working range of the beams. The blue dots show the desired endpoint position for α = ±90°, 0° under 100 N loading. In the isometric view, the undeformed (yellow) and deformed (red) configurations are shown for α = 0° and F = 100 N. The final optimized parameters of all beams are given in Table 1.

Analysis of the objective functions in Table 2 and comparing the errors with Beam 1 as the benchmark shows the effect of tuning different sets of shape-related parameters on achieving desired symmetric behavior in an extended range of motion. Better results were achieved when more parameters were subjected to the optimization process. Furthermore, in Fig. 5 it can be recognized that in comparison with Beam 1, almost all beams exhibited a straighter path of the endpoint subject to a force in the y-direction (α = 0°).

It is fair to note that the anisotropic stiffness in X and Y directions, defined as the range error hasn’t been achieved with circular cross sections. Hence, using non-circular cross sections, or sections with different second moments of inertia about the two main axes, tends to make up for the anisotropy and contributes to favorably agreeable results.

Comparing the results from six optimized beams shows that having more freedom to tune the shape has a lower effect on achieving exceptional results on its own. On the other hand, having more freedom to tune the cross section, effectively enhances the desired behavior. After all, combining these two led to a performance leap to achieve symmetric kinetostatic behavior. It is obvious that a solid conclusion on the exact effect of each optimized parameter needs more investigation of a possibly larger variety of beam parameters.

The logic behind the defined objective function can be used for other applications where anisotropic stiffness and another feature like symmetry are needed. Here, the effect of the two parts of the objective function was chosen to be equal. However, it might be possible to achieve better results for specific requirements with a different weight ratio for f₁ and f₂. Furthermore, each part of the objective function might be improved. For instance, for the symmetry part f₁, it is possible to use more than three paired loadings, which might lead to better symmetry. However, it also has the negative effect of increased computational time.

Our earlier work [29] demonstrated that increasing the number of optimized cross sections does not always lead to better values for objective, while making the problem computationally more expensive. However, the optimal number of control points for the best result may vary in different settings and cases, and this was not the focus of this study. Analysis can be performed to determine this optimal number for each set of requirements.

There are limitations in the self-developed beam model due to its assumptions. Relieving these limitations, suppress all sectional deformations, e.g., warping and in-plane deformations, and using the same parameter optimization approach can lead to more complex deformed shapes and possibly more accurate results. However, advancing the mechanical model has not been the topic of investigation here, and the experiments show that the current model is still valid for the desired loading conditions of this work.

There are several sources of error that might affect the experimental results. Beam fabrication is one of the sources where up to 0.08 mm deviation from the computer-aided design (CAD) model of the beam is observed in printed prototypes. Another source of error is the load angle of the beam, where up to 0.5 deg of error was possible in the fixation of the base rotary stage, which affected the loading angle. The position of the endpoint load was assumed to be fixed at the tip hole of the beam. However, deviation from this location was observed in the different angles of the loading. The correction was applied to the experimental pictures due to camera lens distortion. Still, the image processing and correction might have caused some errors in the measurements. The rotary stage is connected to the base with rigid bars, and processing images in unloaded and loaded states show no measurable deformation of the base. Thus, the beam grounding point can be assumed to be rigid.

Despite the fact that it was anticipated that these two sets of results would differ due to the discussed inaccuracies in the finite element model and potential experimental errors, it is not possible to determine the relative contribution of each to the final deviation of the measurements from the simulation. However, the patterns are identical generally. Consequently, it might be argued that the paper’s message remains legitimate.

It is clear that the resulting beams are not the best possible solution to this problem due to the limited number of available parameters. However, we can state that the presented beam parameter optimization scheme undoubtedly enhanced the desired symmetrical behavior. A similar approach can be applied to an arbitrarily defined kinetostatic behavior for other types of structures with different sets of available parameters.

The combination of parametrization and optimization significantly reduces the input parameters for the whole design process, while using interpolation between those key parameters ensures a smooth transition to avoid discontinuities in the overall structure. Therefore, the proposed method can be assumed as an effective way for the enhancement of other compliant structures’ kinetostatic behavior.

This paper investigates the tuning of different combinations of shape and sectional parameters of asymmetric beams toward enhanced symmetric kinetostatic behavior. An optimization scheme for parametrized beams with an objective function that enhances symmetry in an extensive range of motion was defined. Different sets of parameters representing divergent design conditions were subjected to the optimization. Rather complex shapes were found by this method, and 68% enhancement of the desired behavior was achieved based on experimental results.

The effectiveness of different parts of the method has been validated by comparing the experimental performance of resultant beams from this approach. It is stated that tuning shape and sectional parameters along a slender spatial beam can significantly enhance the desired kinetostatic behavior at the endpoint. The results also suggest solid connections between the availability of shape and sectional parameters for optimization and enhancement of the desired kinetostatic behavior.

The presented method effectively shows its capability to achieve symmetric kinetostatic behavior in the three-dimensional range of motion of the beam’s endpoint. Such design requirements are not easily achievable with existing compliant mechanism design methods. Therefore, this method is an enrichment as a synthesizing method for compliant mechanisms, and it can facilitate the implementation of monolithic spatial compliant elements for other applications.

This work was supported by the Dutch Research Council (NWO) (P16-05 Shell-Skeletons).

There are no conflicts of interest.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.