{'fraction_non_alphanumeric': 0.031007279705156713, 'fraction_numerical': 0.009929639669824252, 'mean_word_length': 4.518910741301059, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Modelling the Nonlinear Response of Fibre-reinforced Bending Fluidic Actuators

Soft actuators are receiving increasing attention from the engineering community, not only in research but even for industrial applications. Among soft actuators, fibre-reinforced Bending Fluidic Actuators (BFAs) became very popular thanks to features such as robustness and easy design and fabrication. However, an accurate modelling of these smart structures, taking into account all the nonlinearities involved, is a challenging task. In this effort, we propose an analytical mechanical model to capture the quasi-static response of fibre-reinforced BFAs. The model is fully 3D and for the first time includes the effect of the pressure on the lateral surface of the chamber as well as the non-constant torque produced by the pressure at the tip. The presented model can be used for design and control, while providing information about the mechanics of these complex actuators.


Introduction
Actuators made by compliant materials are attracting high attention in the recent years [1,2,3,4]. On one hand, soft robotics is proposing a novel approach to build service robots, meaning robots aimed at interacting with humans in their environments [5,6,7]. On the other, soft actuators and sensors find application also in wearable devices thanks to their intrinsic compliance, which makes them interact safely with the human body [4,8,9]. Among all the classes of soft actuators, flexible fluidic actuators (FFA) are one of the most widespread. FFAs are usually composed by a hollow structure made in hyperelastic materials (e.g., silicone rubber) which inflates when compressed fluid is forced into the chamber. We can distinguish two different categories of FFAs based on their actuation modes, the ones that deform axially and the ones that bend.
The most common example in the first category are McKibben actuators [10]. They are composed by a cylindrical ladder, where the compressed fluid (usually air) is inserted, and an external braid, composed by fibers arranged as spirals around the external surface of the cylinder. The consequence of inflation in McKibben actuators is both an axial displacement and a variation of the section. By designing the angle of the fibres it is possible to switch from axial elongation to contraction.
Bending fluidic actuators (BFAs), on the other hand, are slender structures that bend as a consequence of increase in internal pressure. Recently, they started to become widely used, finding applications in legged robots [11], [12], artificial hands [13] and wearable rehabilitation devices [8]. They feature high robustness and easiness in fabrication. An additional advantage is that they fuse structure and actuator in one single element with infinite degrees of freedom but only one degree of actuation (the internal pressure). BFAs can be split again into two groups on the base of the mechanism that transforms the expansion of the chamber into a bending of the whole structure. In the first group the geometry of the chamber is the main element driving the deformation [14], [15], while in the second we have fibre-reinforced bending actuators.
Fibre-reinforced BFAs can be though as McKibbens with the angle of the fibres γ close to π/2, such that the inflation will results in almost pure elongation and no variation in the external perimeter of the cross-section, and with an asymmetric axial constraint ( figure 1 and 2). This constraint, usually a fabric layer glued on one side of the structure, transforms the axial elongation in pure bending with the neutral surface coincident with the fabric layer.
Following the wide adoption of fibre-reinforced BFAs, an accurate mechanical model able to capture their nonlinear response can be very useful to push forward on the adoption of this technology. Modelling soft structures [16] and soft actuators [17], [18] is a challenging research topic due to its multi-disciplinary nature (e.g., fluid-structure interaction) and to the numerous nonlinearities involved (e.g., hyper-elastic materials, large deformations, non-homogeneous properties). Up to know, few models have been proposed in the state of the art. Polygerinos et al. [19] realized a semi-analytical model using a hyper-elastic material model (Neo-hookean). Analytical models have some advantages over FEM models, such as providing insights about the physics of the phenomenon and fast computation. However, the simplifications applied in [19] on the kinematics of the structure make the model practically 1D. The authors did not take into account the effect of the pressure on the lateral surface of the inner chamber, which in facts contributes in the bending due to the incompressibility of the silicone rubber. Finally, they reduced all the effects of the pressure on the tip to a lumped torque constant with the bending angle. The result of all of these simplifications is that their models misses to capture the rich nonlinear behaviour of these actuators.
In this work, we propose an analytical 3D model to reproduce the response of the actuator as a consequence of the increase in internal pressure, when bending in free state. The model is quasi-static but can allows to be expanded to dynamic in further studies. In this model, we take into account the effect of the transverse pressure as well as the variable torque created by the pressure at the tip. Additionally, we analyse all the three principal strains involved in the deformation of the structure and their interactions.
We analysed a slender structure with a semicircular hollow cross-section, as shown in figure 1. The choice of this cross-section is arbitrary (the proposed approach can be applied to different shapes as well) and follows the finding in [19] that this geometry presents a relatively low bending resistance. In order to describe the complex nonlinearities presented by this problem, we decided to use the Eulerian description and apply all the equilibria in the deformed configuration [20]. This approach allowed us to assume a-priori the modes of deformation of the structure and as a consequence apply simplifications that kept the model at an acceptable level of complexity, which enabled the analytical solution of the differential equations involved. The use of numerical computation became necessary only to evaluate algebraic integrals and high degree polynomial equations. As for the description of the kinematics, the use of a toroidal coordinate system allowed us to keep the state of stress principal.
We capture the nonlinear response of the structure and its variation with the design parameter. This result is extremely important because enables the design of the response curve, a clue feature for coupling the actuators with other systems (e.g., the human body) or for exploiting snap-through instabilities [21] or selfstabilized dynamical cycles [12], [22]. Additionally, we analysed the state of stress and strain in the cross section of the structure, which provides important information for design choices about the materials or the dimensions.
Two key features of our model are computational lightness and readability. All the kinematics can be computed once for all and stored in a database, while each step of the statics requires only few numerical integrations, which are computationally efficient. On a standard laptop and using MATLAB for the numerical computations, each step of the statics requires few hundreds of seconds. As for the readability, all of our assumptions are physically motivated and verified a-posteriori. Their validity gives insights into the complex mechanics of these structures and can contribute to the design of novel kind of actuators.
The paper is organized as follows. We first present the design and fabrication of the actuators in section 2. Then section 3 contains the derivation of the mathematical model with all the assumptions, which are verified in the validation in section 5. In the same section we present the results, both in terms of characterization of the response and structural mechanics, and we discuss about their importance. After the circumferential wrapping of the fibres, the fabric layer was incorporated by placing the mock-ups over a thin layer of uncured silicone. (C) Each obtained structure was inserted on a custom 3D printed nozzle that acts as air inlets as well as mechanical clamp. The external diameter of the nozzle is bigger than the internal diameter of the chamber, in order to seal the chamber using the deformation of the silicone itself. With the aim of limiting the movement on a 2D plane without gravity, the experiments where conducted on a lubricated glass plate.
Finally section 6 contains the conclusions.

Design and Fabrication
We produced three different mock-ups of the actuator in order to run the validation tests, whose parameters are contained in table 1. We chose the geometrical parameters to be compatible with the use of the actuators as legs of the soft robot FASTT [12] or as fingers of a glove with rehabilitation purposes [8]. As for the material property µ, we obtained it from a single value declared by the producer of the silicone, rather than by fitting the model over the experimental data. The fabrication of the mock-ups is described in figure 2. The mold to obtain the main body of the actuator has a semi-cylindrical shape and an insert of the same shape to create the internal chamber (figure 2 A). Both the main mold and the insert were 3D printed in ABS plastic. We used as silicone rubber the Smooth-on Ecoflex 0030 (tensile modulus at 100% strain E 100% = 69 kPa, ultimate strength σ Rs = 1.38 MPa, elongation at break 900%). We obtained the value of µ in table 1 simply as µ = (4/3) * E 100% . As for the fibre, we used a wire featuring high tensile stiffness and strength (declared ultimate strength σ Rw = 5.8 GPa), but negligible bending stiffness (diameter d = 0.10 mm) to constrain the radial expansion. A fabric with squared stitches made in high-density polyethylene (HDPE), with a tensile modulus E f = 1.5 MPa was used to constraint the axial elongation (figure 2 B and C).

Mathematical Model
The goal of this model is to describe both the response of the actuator and its mechanics in terms of stresses and strains, capturing all the nonlinearities deriving from: (1) the hyper-elastic material, (2) the large deformations, (3) the presence of the fibre reinforcements. We wanted to keep the model easy to read and computationally light. The readability allows the manipulation of the model to adapt it to different structures and provides hints about the physics of the actuator, which can be used to improve its design. As for the computational cost, a light model is suitable for control purposes, a highly required feature in robotics applications. For all the above reasons, we decided to make simplifying assumptions based on hypotheses on the mechanics of the actuator, which allowed us to obtain a semi-analytical model, where we analytically solve all the differential equations describing the kinematics and the equilibrium of the structure.
The first assumption that we make is to consider the deformation of the structure to follow a toroidal symmetry, as described in figure 3. In particular, we divide the cross section in three main parts (figure 1): (1) circular, (2) rectangular, (3) corners. We assume the circular part to deform as a portion of torus, while the rectangular part and the boundaries as a portion of cylinder. This assumption implies to have the curvature of the actuator to be constant in each point along the axis and equal to χ = 1/R, corresponding to the curvature of the surface coincident with the inextensible bottom layer (figure 1).
The model follows the Eulerian description, so we assume a-priori the deformed configuration and we write all the equilibria in this configuration. We then use these equilibria to find the value of the internal pressure that can satisfy them. As far as the nomenclature is involved, we used a star . * to indicate quantities in the undeformed configuration that change in the deformed one.
All the assumptions together with their motivations are contained in the remaining of this section, which is organized as follows: we treat individually the three portions and for each of them we first describe and solve the kinematics, then we include the material model and finally we apply the equilibria and solve the equations. η r θ ϕ z Figure 3. Description of the structure deformed following a toroidal symmetry, with the associated toroidal coordinate system r, θ, φ.

Semicircular part
3.1.1. Kinematics We assume the circular portion of the structure (1 in figure 1) to deform into a section of torus defined as (figure 3) In order to define the kinematics of this structure, we reduce the problem to just a half of the section due to symmetry (0 ≤ θ ≤ π/2) and we introduce nondimensional parameters to describe the geometry: Because the curvature χ stays constant along the axis, all the quantities are the same in each cross section ( ∂ ∂ϕ = 0); on the other hand, the particles deform along ϕ due to the bending. We can describe this deformation by means of λ ϕ Now, in order to describe the kinematics along the r and θ directions, we need to introduce two functions g and f , which map the values of the radius r and angle θ for a material particle in the deformed configuration back to the values r * and θ * in the undeformed configuration. We will define these functions referred to the non-dimensional parameters Now we can introduce the strains along r and θ, relating them to these functions Let know apply the incompressibility constraint, obtaining the linear PDE that describes the kinematics of the system Considering that f (ρ, θ) has to fit two boundary conditions on θ: we assume the displacement of the particles along θ to be negligible and as a consequence we set f (ρ, θ) = θ and thus ∂f (ρ, θ)/∂θ = 1. Substituting in 5 we obtain which we can integrate respect to ρ to obtain By applying the boundary condition on ρ (g(1, θ) = 1, meaning that the external surface remains undeformed) we can find c(θ) Finally we find an algebraic expression for g(ρ, θ) where we excluded the negative case because g(ρ, θ) ≥ 0 by definition. From (10) we can see that when the actuator is at rest (Γ = χR 2 = 0), then g(ρ, θ) = ρ , ∀θ ∈ [0, π]. As far as the conditions of existence for (10) are concerned, we need Physically, this condition fixes the minimum limit for the deformed internal radius R 1 , which is the minimum value of ρ, in order to have an undeformed radius R * 1 ≥ 0. This means that we can make the condition strict.
Finally, using (8) the value of the lower boundary for ρ, i.e., k = R 1 /R 2 can be computed by solving: From (12) we can see that: (1) we can compute k knowing the design parameters and the value of χ; (2) k depends on θ, which means that in the deformed configuration, in order to respect the conservation of volume, the arcs formed by the particles in the reference configuration will deform and not stay circular any more; (3) if ρ ∈ [k, 1], the condition of existence (11) will be automatically satisfied. The kinematics of the semi-circular part is now solved, so for each assumed value of χ we can compute all the strains (4) at each point (ρ, θ, φ).

Equilibrium
The equilibrium in the deformed configuration has to be expressed in toroidal coordinates, coherently with the description of the kinematics. In this configuration, because the only external action was provided by the internal pressure, we assumed the state of stress in the deformed configuration to be principal, obtaining We can now apply the Eulerian field equation in static equilibrium and with no body force ∇ · T = 0 [20], which expressed in toroidal coordinates from [23] gives It is interesting to observe that the equilibrium in the radial direction (14) contains also the stressesσ φ and σ θ , as a consequence of the bending of the longitudinal axis from the undeformed cylindrical configuration to the deformed toroidal one. This fact confirms the presence of a coupling between the three principal stresses.
Due to the assumptions that we made on the kinematics, the circumferential (15) and the axial equilibria (16) become redundant, so we will apply only the equilibrium (14).

Material Model
As far as the material model is concerned, we wanted to find a trade-off between accurateness in the description of the behaviour and compactness, to allow easy manipulation in an analytical model. We discarded a linear material model because one of the key points of this work is to capture the nonlinear response of the actuators. We considered then the Neo-hookean material model [20] and the Gent model [24]. The first one has a relatively compact formulation and can replicate accurately the initial response of hyper-elastic rubbers, while it totally fails in capturing the stiffening that happens when the molecules become almost totally stretched. The Gent model, on the other hand, can accurately capture all the nonlinear response of rubbers but it has a slightly more complicated formulation. By comparing the two models, we could observe that the discrepancy in the behaviour happens for λ > 3. From a first estimation we expected the strains to be λ ≤ 3 for our actuators in the chosen ranges of design parameters, so we decided to use the Neo-Hookean model and then to verify aposteriori the goodness of our choice. The constitutive equations for the incompressible Neo-Hookean materials are then where Q(ρ, θ, χ, p) =Q/µ is a nondimensional lagrangian parameter.

Solution
We can rewrite equation (17) as with a 1 (ρ, θ) = 1 + 2Γρ sin θ/(ρ(1 + Γρ sin θ)) (20a) a 2 (ρ, θ) = Γ sin θ/(1 + Γρ sin θ) Substituting (18) into (19), we obtain Equation (21) is a linear first order differential equation, which solution is We need now a boundary condition on the radial stress σ ρ , which we can apply at the inner chamber where it has to equal the internal pressure Substituting equation (18) into the boundary condition (25), we get The problem for the semi-circular part of the cross section is now closed: knowing χ and the strains obtained from (4), we can compute all the stresses (18) for a given value of the internal pressure p.

Kinematics
The rectangular part represents the wall between the chamber and the bottom layer (2 in figure 1). We assume that this portion of the section, whose domain is defined in (27), deforms as a consequence of the internal pressure, while the corners at the boundaries (3 in figure 1) remain unchanged. We will treat the kinematics of the corners in section 3.3. The rectangular portion when deformed becomes a sector of a hollow cylinder defined as (figure 3) Non-dimensionalizing the variables and considering only half a cylinder due to symmetry we get In this domain, we can define the strains. The longitudinal strain is For the transverse strain, due to the constraint represented by the external fibres on both sides of the structure, we assume there is no movement of particles in the ζ direction and so λ ζ = 1. Now, we can obtain λ η from the conservation of volume λ ϕ λ ζ λ η = 1 ⇒ λ η = 1/λ ϕ = 1/(1 + Γη).
In order to find the deformed thickness β we can integrate λ η (31) Solving we get (32) β = ( 1 + 2Γβ * − 1)/Γ were we select only the positive solution because we want β ≥ 0. With this result we completed the kinematic analysis and we can obtain all the strains in the rectangular part (30) for each given value of χ.

Equilibrium The only non-trivial equilibrium equation expressed in cylindrical coordinates is
which transformed in non-dimensional form becomes where G(η, χ, p) =G/µ is a nondimensional lagrangian parameter.

Corners
For a matter of simplicity, we restrict the kinematics of the corners (3 in figure 1). In particular, we assume that for the corners the conservation of the volume does not hold. This statement may seem physically not plausible, but we should consider that both the circular part and the rectangular part will in fact move some of their particles into the corners. But, if we followed these particles, the kinematic description for both the circular part and the rectangular part would be much more complicated than what we obtained in the previous sections. In order to keep an acceptable level of accuracy and at the same time to obtain an elegant mathematical model, we then decided to simplify the description of the corners.
With these assumptions, each of the corners deforms as a portion of hollow cylinder. We consider only one for symmetry, which description is Non-dimensionalizing the variables we get Following the previous assumptions we obtain λ ζ = 1, λ η = 1. The only non-zero deformation is then λ ϕ = 1 + Γη.

Rotational equilibrium
In a first model for BFAs, Polygerinos et al. [19] reduced all the effects of the internal pressure p to a concentrated torque M t acting at the tip of the actuator, generated as the integration of the infinitesimal torques on the transverse section at the tip respect to the inextensible bottom layer: where in our case R 1 = R 1 (θ). From this equation, we can notice that M t depends linearly on p but does not depend on the curvature χ. While the torque M t can be considered the main action that produces the bending when the actuator is at rest (χ = 0), it fails to capture its behaviour when the deformation increases. In facts, from experimental results we obtained that these actuators show a highly nonlinear response for 0 ≤ χL ≤ 3π/2 (the typical range of application) as shown in figure 5. One essential element lacking to reproduce this trend is the pressure force at the tipF t = 0.5pπR 2 1 , which generates a moment M χ respect to each section ϕ. M χ , non-zero for χ > 0 and depending on χ, respect to a generic section ϕ equals (see figure 1) Summarizing, the pressure at the tip of the actuator produces a force F t , directed as the tangent at the tip and applied at the centroid of the transverse section of the chamber; however, because the neutral surface is forced to coincide with the inextensible layer, the transport of F t to the neutral surface produces also the moment M t . Now, the non-dimensional internal elastic torque, constant in each of the sections ϕ, is where apexes c, r, and b refer to circular part, rectangular part and boundaries (corners) respectively. Let now non-dimensionalize and compute the external torques. For the transport torqueM t , first of all we need to add the contribution of the area created by the compression of the rectangular part and then we can non-dimensionalize dividing by µR 3 2 and integrate along ρ where the last integral has to be computed numerically due to the presence of k(θ). Non-dimensionalizing the moment of the force F t (49), we have instead For the scope of our work we assumed the actuator to deform with constant curvature, so all the torques will be the same in each section, but in reality M χϕ depends on ϕ. For this reason, we decided to apply all the equilibria in the middle section, in which ϕ = 0.5, so where again the last integral with k(θ) has to be computed numerically. At the end, our equilibrium equation becomes

Solution
We have now all the equations to describe the kinematics and statics of our actuator. Interestingly, we could solve analytically all the differential equations involved in the model, while we need to use numerical computation to evaluate some algebraic integration. The procedure that we decided to follow to solve the equations is the following. We set a value of the curvature χ = Φ/L, being Φ the angle between the tangent at the base and the tangent at the tip of the actuator ( figure 1). The range of Φ that spans all the physically possible configurations of the actuator when bending in free state is [0, 3π/2]. For each χ, we can compute all the kinematic quantities including strains and displacements. In order to compute the quantities that vary on the plane of the cross section, we created three equally spaced grids of points for ρ and θ in the circular section and η and ζ in the rectangular section and in the corners ( figure 6). Next, we can obtain the solution of the statics by starting with a guessed value of the internal pressure, called p g . Using this value and the already obtained kinematics, we can compute all the stresses on the section using the same grids defined for the kinematics and finally we can compute the error function for the equilibrium, definite as F (χ, p g ) = M σ + M t + M χ . Using an iterative algorithm we can now solve the problem of finding the value of p that guarantees the equilibrium, for each χ, by setting The model and its solution are now closed. In the next sections we will describe the experimental set-up (section 4) used to validate the model (section 5.2) and the results on the mechanics of the structure (section 5.3) and on its characterization as a soft actuator (section 5.4).

Experimental set-up
The experimental set-up (Fig. 6) was composed by: (1) one proportional electronic valve (Camozzi K8P), Input: 0 − 10 V, Output: 0 − 3 bar; (2) an Arduino MEGA microcontroller; (3) a camera for data acquisition; (4) the mock-up. Each experiment was conducted starting from rest position of the actuator and increasing the pressure up to the desired value through the electronic signal provided by the Arduino microcontroller to the electronic valve. The actuator bended as a consequence of the increase in internal pressure and its equilibrium configuration was recorded using the camera. From the pictures we computed the bending angle using Kinovea software for image analysis (figure 4), with a resolution of 1 • . Each experiment was repeated 4 times and we computed the average and the standard deviation.
A big advantage of this model on the side of computational time is that all the kinematics can be solved once for all for each values of χ and the results can be stored in a database. As for the statics, using 610 points on the grids of the cross-section, MATLAB and a simple algorithm (not optimized for fast computation) running on a 2.9 GHz Intel Core i5 computer, each step of the statics is obtained in ≈ 0.05 s.

Validation
In order to validate the model we compared its results with the experimental ones. As for the model, we span the Φ angle in [0, 3 2 π] and for each value we obtained the pressure p that guarantees the equilibrium of the structure. The experimental results were instead obtained by driving the actuator with a pressure p in [0, 0.35] bar and acquiring the correspondent angle Φ at the equilibrium. The validation of the model on the three mock-ups is shown in figure 5. It is possible to observe that the behaviour of the actuator changes considerably between the three mock-ups and that the model is able to capture the response in all of the cases. It should be noticed that the validation is successful using the nominal parameters for the geometry and the material, without any fitting procedure, confirming the physical plausibility of our assumptions.
We wanted also to explore the importance of the pressure on the lateral surface of the chamber respect to the one on the tip. We then artificially set to zero the pressure on the lateral surface and simulated the model in the same interval. The results can be found in comparison with the full model and the experiments in figure 5 under the name "model pl OFF".

Strains and stresses
In this section we analyzed the state of deformation and stress on the cross sections of the structure. We used the same parameters of the mock-up number 2 (table 1) for a bending angle Φ = 3/2π. In figure 6 we can visualize the deformed cross section, with the equallyspaced grid used for the numerical computation, overlapped over the undeformed cross section. The deformation of the chamber is clearly macroscopic and follows the assumptions made about the kinematics in sections 3.1.1, 3.2.1 and 3.3. As assumed, we can observe a sliding between the rectangular part, the corners and the circular part, which we used as a simplification of the behaviour of the boundaries that  connect the three different portions of the cross section. We computed the maps of strains and stresses on the cross section of the structure. The results of the stresses and strains along the rectangular coordinates η and ζ are shown in figure 7.
As we can notice, due to the assumptions there are little discontinuities in some of the values of stresses and strains, mainly localized around the corners where the three sections touch each other. As for the strains, it is interesting to notice that the maximum strains are all ε ≤ 1.2, so λ ≤ 2.2, thus confirming the formal equivalence in the use of the Neo-Hookean model respect to the Gent model, as it was assumed in section 3.1.3. We can notice that the strain in longitudinal direction ε ϕ , which is the one driving all the others, is homogeneous over all the cross-section and increases continuously with η. As expected the For a better understanding of the stresses and strains behaviours on the circular part of the cross section, we also plotted the ρ and θ components in figure 8. We can observe that, coherently with the assumptions of section 3.1.1, ε θ and ε ρ depend on both ρ and θ. We have a state of compression along ρ while a state of tension along θ. In both cases the maximum values are localized at the internal radius ρ = k(θ).
The results on the response of the actuator allow the functional design of the actuator. In facts, from the results in figure 5 we can choose the value of static and dynamic pressures (p 0 and ∆p) in order to obtain a certain target oscillation Φ 0 + ∆Φ. Additionally, by knowing the level of stress and strain we can design the structure in terms of materials and dimensions, fixed the value of the pressure and we can predict where the structure will eventually fail.

Characterization of the actuator
We computed the angle-pressure characteristic for different values of the k * parameter in order to characterize the actuator. The parameters used for the computation correspond to the mock-up number 1 (table 1) but with R * 1 changing accordingly to k * , while the results are shown in figure 9. As we can observe, the response of the actuator changes gradually by varying k * , so it is possible to tune the response of the actuator by choosing the design parameter k * .
We then chose two values of the ratio k * and analyzed the characterization for different values of q * = L * /R 2 . The results are shown in figure 10. The parameters are the same used before (table 1, mock-up  1), but with R * 1 and L * varying according to k * and q * respectively. In this case, as expected, q * mainly influences the reactivity of the actuator, which results more steep for more slender structures (q * high) and more gradual for lower values of q * .
Apart from the mechanical design of the actuator, the last two results obtained in this section are extremely useful for the functional design of the response of the actuator. In facts, being able to tune the nonlinear response of soft actuator is a highly beneficial feature, allowing the design of systems presenting snap-through instabilities to amplify the response [21] or self-stabilized limit cycle oscillators [27,12]. On the other hand, by choosing a high value for k * , which corresponds to a relatively thin wall thickness, the response can be close to linear.

Conclusions
In this work, we presented a quasi-static, 3D mechanical model for fibre-reinforced BFAs. The model includes the effect of the pressure on the lateral surface of internal chamber of the actuator and the non-constant torque induced by the pressure at the tip. All of the three principal strains and their interaction are taken into account.
In order to keep the complexity of the model at an acceptable level, we used the Eulerian description and applied all the equilibria in the a-priori defined deformed configuration. The resulted formulation is readable and computationally efficient: it is possible to compute the kinematics once for all and store the results in a database; each step of the statics requires only numerical integrations of algebraic expressions.
We experimentally validated the model and then used it to reproduce the nonlinear response of these actuators and its variation with the design parameters. This mechanical model can be useful in several ways: (1) it is a powerful tool for designing the actuator response for a target application; (2) can be used for the mechanical design of the structure (e.g., materials, dimensions); (3) by capturing the intrinsic mechanics of these systems, provides insights for the design of novel types of soft actuators.