{'fraction_non_alphanumeric': 0.03914698140727949, 'fraction_numerical': 0.012691334015122303, 'mean_word_length': 4.340974212034384, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Mechanics-based Analysis on Flagellated Robots

We explore the locomotion of soft robots in granular medium (GM) resulting from the elastic deformation of slender rods. A low-cost, rapidly fabricable robot inspired by the physiological structure of bacteria is presented. It consists of a rigid head, with a motor and batteries embedded, and multiple elastic rods (our model for flagella) to investigate locomotion in GM. The elastic flagella are rotated at one end by the motor, and they deform due to the drag from GM, propelling the robot. The external drag is determined by the flagellar shape, while the latter changes due to the competition between external loading and elastic forces. In this coupled fluid-structure interaction problem, we observe that increasing the number of flagella can decrease or increase the propulsive speed of the robot, depending on the physical parameters of the system. This nonlinearity in the functional relation between propulsion and the parameters of this simple robot motivates us to fundamentally analyze its mechanics using theory, numerical simulation, and experiments. We present a simple Euler-Bernoulli beam theory-based analytical framework that is capable of qualitatively capturing both cases. Theoretical prediction quantitatively matches experiments when the flagellar deformation is small. To account for the geometrically nonlinear deformation often encountered in soft robots and microbes, we implement a simulation framework that incorporates discrete differential geometry-based simulations of elastic rods, a resistive force theory-based model for drag, and a modified Stokes law for the hydrodynamics of the robot head. Comparison with experimental data indicates that the simulations can quantitatively predict robotic motion. Overall, the theoretical and numerical tools presented in this paper can shed light on the design and control of this class of articulated robots in granular or fluid media.


I. INTRODUCTION
Apart from animals capable of crawling, digging [1], slithering [2], swimming [3], and gliding [4,5] in environments, bacteria, one major group of microorganisms, also inspire the development of novel and efficient robots. 90% of marine bacteria [6] achieve efficient locomotion in a fluid through the rotation of a flagellum (a slender filament), exploiting the anisotropic drag to produce thrust and violating the constraints of the scallop theorem [7]. Flagellar locomotion results from a non-trivial coupling between the geometrically nonlinear deformation in the flagellum and mechanics of the surrounding medium, posing a challenging fluid-structure interaction (FSI) problem.
Theoretical study about flagellar propulsion dates back to 1955 when Taylor [8] first analyzed the swimming of microscopic organisms. Over the last two decades, numerous studies have studied flagellar propulsion in low Reynolds fluids through experiments [9][10][11][12][13][14], computation [15][16][17], and theory [18][19][20][21]. Recent investigations [22,23] have modelled the flagellum as a Kirchhoff elastic rod [24], and coupled to the fluid with hydrodynamic forces [25]. Jawed investigated the dynamics of a helical elastic flagellum rotating in a viscous fluid [26] and near a rigid boundary [27]. However, the * Correspondence should be addressed to khalidjm@seas.ucla.edu role of the head and the flow generated by its motion and coupling with the flagellum-induced flow are ignored. Huang performed simulations that utilized flagellar buckling to change moving direction on a robot composed of a mass-point head and uni-flagellum [28]. Our previous work [23] established an untethered articulated robot that was composed of a rigid head and multiple soft flagella and used discrete differential geometry to simulate the flagella and resistive force theory [25] (RFT) to model the interaction force between the fluid and flagella. Simulation and experimental results agreed well quantitatively. Our robot also demonstrated the same behavior as bacteria, ii.e., head and flagella rotating in the opposite directions [29] and circling when near the air-fluid boundary.
RFT, initially used for analyzing the movements of microscopic organisms in viscous fluids, is proved to apply to animal and robot locomotion on and within granular media (GM) [30]. This builds an intimate link between the microscopic bacterial world and meter-sized animals in sand. The most recent studies [31,33], on rotational intruders moving through GM conclude on the empirical feasibility of applying modified and entirely empirical RFT, e.g. granular RFT [33]. These experimental studies demonstrate that RFT is simple yet very effective in granular materials considering the complicated constitutive features of granular matter such as nonlinearity and nonlocality. However, the intruders in these works are rigid, inspiring our further study on the geometrical arXiv:2112.13927v1 [cs.RO] 27 Dec 2021 design and efficiency study on robots with soft flagella in GM. Previously, our other work modeled the locomotion of articulated soft robots in GM [22]. Leading design and control parameters of our untether robots are the rotational speed of the embedded motor, the number of flagella, geometrical parameters of flagella such as radius and length. Some counter-intuitive observations include the inverse relationship between the speed and the number of flagella of the robot in the representative setup (as shown in Fig.1(a1-a5)). This verifies that the flagellar locomotion is the result of a complicated coupling between granular mechanics and deformable bodies, a complex FSI problem.
In summary, while there exist works [13,23,28] that explore flagellar locomotion considering the effect of the head, systematic simulations are required for more quantitative predictions, e.g., under what design and control parameters the maximal efficiency of actuation is achieved. Furthermore, although a comprehensive simulation framework is offered in our earlier work [22], elaborate equations suffocate and make it difficult for designers to understand and predict the movement of the system directly. As a result, a reduced model capable of qualitatively capturing the relationship between the performance of flagellar locomotion in GM versus the design space is required, hence avoiding the cumbersome trial and error design process.
Here, we employ a structural robot design similar to that used in our previous work [22]: a palm-sized untethered robot composed of n ≥ 2 naturally straight elastic rods and a rigid head housing a geared dc motor and batteries. As illustrated in Fig. 1, the rotation of flagella induces the drag force from the surrounding GM because of their flexibility. Flagellar deformation provides a net propulsive force forward, which is zero in the case of straight rigid flagella. Meanwhile, we implement two numerical frameworks, both of which employ Resistive Force Theory (RFT) to model the drag on a collection of multiple elastic rods and modified Stoke's law to calculate the drag force and moment on the head. In contrast, one uses Euler-Bernoulli beam theory, and the other applies discrete differential geometry (DDG) to simulate the flagellated robot structure. When the motor rotates slowly, e.g., ω T ≤ 10 rpm, the deformation of the robot flagella is linear, and the resulting propulsive force calculated using DDG simulation is nearly identical to that obtained from Euler-Bernoulli beam theory. However, the beam theory fails to precisely capture the nonlinear deformation of the robot flagella at high motor speeds, e.g., ω T = 250 rpm, DDG-based framework is still capable of accurately representing the performance of the robot in experiments. Due to the rod-based kinematic representation of the robot, this computational tool is used to simulate flagellated robots even faster than real-time. Additionally, our prior study uncovered a counterintuitive phenomenon, the speed of the robot falls as the number of flagella rises (design 1), which is captured by both the beam theory and DDG-based frameworks. However, the beam theory-based framework predicts the existence of design 2, in which the speed of the robot increases as the number of flagella grows. This prediction is successfully confirmed experimentally after we modify the robot design and radius of the GM. In summary, the simple overall Euler-Bernoulli beam theory-based framework can qualitatively capture both designs of the complicated locomotion in GM, the intricate interaction between the robot and GM. In contrast, the DDG-based framework is capable of quantitatively simulating the complex locomotion in GM.

II. PROBLEM STATEMENT
In this work, we analyze flagellated locomotion in GM, a typical FSI problem, using a mechanics-based approach. As a result, both the experimental and theoretical aspects of the problem involve two primary components, the GM (fluid part) and the robot (structure part). A GM must be selected to complement the architecture of the untethered flagellated robot. Internal friction between granules prevents the robot head from fluidizing the granules in front and propelling forward, whereas insufficient friction results in ineffective flagellar propulsion. Internal friction during locomotion is correlated with the phase transition of GM from solid to fluid, which is governed by temperature and volume fraction (VF) [44] -the ratio of solid to occupied volume. The less effort required to modify the VF of GM, the better for investigating the locomotion in GM. Apart from the volume fraction of the GM, other physically controllable parameters include the shape and material friction of the head, as well as the number, radius, length, and stiffness of the flagella and the rotational speed of the motor. Among them, the number, n, and the rotational speed of the motor, ω T , are the most precisely trackable with the least effort owing to the design simplicity of the robot and experiments. Their effect on the translational speed of the robot, v, is explored and will be displayed in Sections V A and V C 1.
The entire robot structure consists mainly of active multiple soft elastic flagella and a rigid passively actuated head, all of which rotate along the long axis of the robot. The case where the flagella are fixed at a 3Dprinted plate with uniform external drag forces applied is analogous to the case of a cantilever beam with a uniformly distributed load. In what follows, first, we outline the Euler-Bernoulli beam theory-based mechanics analysis of the untethered robot, including (1) the external loading from the GM onto flexible flagella and (2) the drag force on the head. Beam theory qualitatively captures the two cases in which increasing n can either accelerate or retard the speed, v. Whichever happens, in reality, is determined by the intricate balance of the competition between the external loading and elastic forces. According to experimental evidence, the result is closely related to the robot head and motor speed design, ω T .
Second, we introduce a numerical model of the robot in which the robot structure is represented by a network of Kirchhoff's rods[]. Finally, we present experiments conducted to quantify the propulsive speed of this class of flagellated robots in Section IV.

A. Beam theory based analysis of propulsion
Benefiting from the elasticity property, the soft flagella of the robot bend out of the long axis, and the drag force generated by the surrounding granules produces a component in the -y direction as displayed in Fig. 2, propelling the robot forward.

Clamped flagellum fixed in space
For illustration, one flagellum of the robot in a discrete setting moving in the GM is displayed in Fig. 2. Hereafter, unless otherwise stated, all the parameters, e.g. forces, are associated with one flagellum. Consider a node x(s) along the flagellum, where s ∈ [0, L] is the arc-length parameter and L is the length of the flagellum. At steady state, the magnitude of the velocity vector at where v is along −y axis, ω t stands for the rotational speed of flagella and R d is the radius of 3D printed plate holding soft flagella. The velocity, v at this point can be decomposed into two parts: the parallel term v t = (v · t)t = v t t and the perpendicular term v where v t = −v sin θ t and t represents the tangent at the point, whose direction is shown in Fig. 2. Inspired by Maladen et al. [30] and Texier et al. [31], the drag force exerted by the GM on the robot is modeled using v h x s=L s=0 Schematic showing the velocities and induced resistive forces on one flagellum rotating at a constant rotational speed ωt on the robot moving forward at a velocity of v h . The coordinate system (O − x − y − z) is fixed at the center of the 3D-printed plate; the flagellum starts at s = 0, propagates along the x-direction until the end (s = L), and is divided into a series of nodes. Each node (numbered based on differential geometry) is characterized by its tangential direction t and its velocity v direction; k-th node (with coordinate x(s)) experiences a force F ext k that propels the robot forward (along −y-direction).
generalized Coulomb's friction law in which the friction force experienced by a slender cylinder is broken down into components normal and tangential to the segment axis, with two corresponding force coefficients. The drag is shown to be dependent on the shape of the object, e.g. the diameter [32] . Additionally, one work [46] considers the effect of the geometry of the soft filament while modeling the local drag. As a result, the tangential and perpendicular force constants (force per unit length) from the medium that resist v t and v p are where the drag coefficients, one along the tangential direction and one along the perpendicular direction are where µ is the coefficient constant that quantifies the robot-granule friction as a function of the granule size and the inertia and surface friction of the head, and L and r 0 are the length and radius of each flagellum, respectively. The drag force constant of external forces at node x(s) is Then, the component of F ext along the x-axis (i.e. antiparallel to velocity) is where sin θ = dw and cos θ = dy The transverse displacement of a flagellum, which is modeled as a beam is w(s). Assuming a small deflection, w(s) ≈ w(y) where w(y) represents the flagellar deflection along the y-axis. We ignore the effects of head translation in this section; they will be discussed in the following section. The viscous force can be split into y-and xcomponents; the former produces propulsive force, while the latter produces a torque that rotates the head. Given our definition of the direction of y-axis, the force constant of propulsive force is as follows along −y: Without ignoring any higher order terms, i.e. sin 2 θ, the x-axis force constant is formulated: where EI is the bending stiffness of the beam (i.e. the flagellum).

Clamped flagellum moving at constant speed
During the robotic movement, in addition to the flagella, its head is also rotating and translating. Assuming the robot moves along the y-axis at speed v h , we rewrite its velocity as wherex,ŷ are the unit vectors in the direction of the x-axis and y-axis, respectively. The force constant of viscous drag force is: and its component along y-axis is Adding q h to Eq. 8, we get the final viscous force We integrate this across the length of the flagella to compute the total propulsive force. Its component along the x-axis is Adding this solution to Eq. 9, we get the following expression for the force constant along x-axis: We used the bvpfcn function in MATLAB to solve this nonlinear fourth order differential equations with boundary conditions where the prime ( ) indicates differentiation with respect to y. Since the head-induced effect is taken into account and no approximations about the flagellar deflection are made during the calculation, this design is referred to as nonlinear beam (NLB) regime. However, by making some simplifying assumptions, we can achieve a more straightforward, closed-form solution. If we assume that the flagellum is slightly deflected to be treated as an Euler-Bernoulli beam, then sin 2 θ ≈ 0 and cos 2 θ ≈ 1. This enables us to simplify the solution and achieve linear approximation, referred to as linear beam (LB) regime. The curvature in Euler-Bernoulli beam theory is approximated as d 2 w dy 2 and dw dy 2 terms are assumed to be negligible compared with 1. The boundary conditions on the beam are Eq. 17. The deflection of a LB is, therefore, calculated as follows and p t here is obtained as follows by omitting the sin 2 θ term in Eq. 9: Hence, the simplified total force along x-axis is The linear propulsive force is likewise simplified as and the total propulsive force is Note that we can neglect the head velocity in these equations because it has been experimentally proved to be small. The deflection of LB, NLB, and NLB w/o head (nonlinear beam without head) is compared in Fig 7 and Section V D. In addition, the head experiences a resistive force due to its translation and rotation. The viscous drags due to its translation and rotation are respectively, where 6πµR h v h is the drag force and 8πµR 3 h ω h is the torque on a perfectly spherical object in a low Reynolds number fluid according to Stoke's law and v h =ẋ 1 is the velocity of the head, and ω h =θ 0 is the rotational speed of the head. The coefficients C 1 and C 2 are included because our robot head is not a sphere; given the difficulty of measuring C 1 , C 2 , and µ during experiments, they are used as fitting parameters (see Section V B) to match the experimental and simulation results in Section V. Furthermore, this flagellated robot system is balanced in terms of both force and torque. As a result, the propulsive force generated by the flagella should be equal to the drag force on the head, F p = F p , and the torque actuating the rotational movement of n flagella should also be equal to the moment on the head. Using F p = y=L y=0 q(y)dy, or the corresponding linear approximation in Eq. 22, we can compute the velocity: where R h is the radius of the robot head. Next, we calculate the rotation speed of the head, ω h , based on the torque balance. To compute this, we need the torque generated by F x , which is given by where F x = y=L y=0 p t (y)dy, or with the linear approximation from Eq. 20. Given that the torque generated by the flagella (in Eq. 26) must be equal and in the opposite direction of the torque on the head in Eq. 24, we can compute the rotational speed of the head: Plug Eq. 6 into Eq. 9, we obtain p t as a function of only y and ω t , denoted as p t = f (y, ω t ) and thus F x = f (ω t ). Then, from Eq. 27, we find that ω h is a function of ω t whereas the remaining parameters are constants. In addition, one of control parameters in experiments is the rotational speed of the motor ω T : Discrete schematic diagram of a robot with two flagella. (a) Geometry of the robot in its undeformed state. Here, L h = 2R h is the diameter of the robot head, Lp is the diameter of the 3D-printed plate bridging the head and flagella, and L is the length of each flagellum. Dashed lines represent the rigid structure, while solid lines represent the elastic structure. The node x1 specifies the center of the head, i.e. its location. (b) A general close-up view of three adjacent nodes, x k−1 , x k , and x k+1 , and two edges, e k−1 = x k − x k−1 and e k = x k+1 − x k . The turning angle, ψ k , between the two edges results in bending energy and the rotation of the material frame from one edge to the next results in twisting energy. The reference frame on e k is d k The twist angle from the first material director associated with t k−1 to the next material director associated with t k is θ k . (c) A close-up of the "T"-shape joint node x2 that connects the head to the flagella. d.B.E. indicates discrete bending and twisting energy. Only the joint node, x2, is connected to more than two nodes.
Note that ω h and ω t here are scalars, but the rotational directions of the head and flagella of the robot are opposite. In conclusion, we have two equations for two unknown variables ω t and ω h , allowing us to solve the propulsion of the entire system. From these two variables, all other parameters, such as v h and F p (by combining Eqs. 22 and 25) can be calculated. We emphasize that the overall framework includes not only the flagellar propulsion based on beam theory but also the head motion based on Stokes'law.

B. Numerical method of flagella locomotion
In this section, we present the numerical simulation framework in which the robot structure is modeled as a network of connected Kirchhoff elastic rods. The elastic energies of this structure are thereby the linear sum of bending, twisting, and stretching energies, the negative gradient of which is the elastic forces on each degree of freedom (DOF). The implicit Euler method is used to solve equations of motion (EOM) in which the external drag forces from the GM on each DOF are modeled through RFT. When ω T is small, locomotion parameters obtained from this simulation, e.g., velocity, propulsive force, and deflection, are compared to those obtained from analytical equations based on Euler-Bernoulli beam theory described in Section III A. Furthermore, the simulation results of DDG are shown to match quantitatively with experimental data, validating the method itself.

Discrete Differential Geometry (DDG)
Fig 3 is a discrete representation of a robot with a rigid head and 3D-printed plate and n straight elastic flagella (n = 2 in the figure). The numerical values for physical parameters will be given in Section III C. Along the schematic discretization, there are a collection of N nodes (circles in Fig. 3 , called edges (the lines between two adjacent nodes). Hereafter, we use subscripts for node-based quantities and superscripts for edge-based quantities. As shown in Fig. 3(b), the head is discretized into x 0 , x 1 , x 2 , with node x 2 denoting the one at the head-flagella junction. It is a unique node that is connected to n + 1 additional nodes. All other nodes are connected to two nodes, or to a single node in case of the last node at the end of each flagellum, denoted by an open circle.
Rather than being just a single vector, each edge (see Fig. 3) is also equipped with two corresponding sets of orthonormal reference frames to track its rotation, a reference frame d k 1 , d k 2 , t k and a material frame m k 1 , m k 2 , t k ; both of them share the tangent t k = e k / e k as one of the directors ( · represents the Euclidean norm of a vector). The reference frame serves as a frame initialized at time t = 0 and updated at each subsequent time step via time-parallel transport, as illustrated in Fig. 3(c), the material frame can be calculated in terms of a scalar twist angle, θ k . The detailed transformation expression between the reference and material frames can be found in Du et al. [22] and Jawed et al. [34]. For this class of flagellated robots, the DOF vector is constituted of node positions and twist angles formulated as follows (29) where N e is the total number of edges in the entire robot, which in this work is N − 1 (see Fig. 3(a)), and the superscript T stands for transpose. q has a dimension of ndof = 3×N +N e . Since q completely defines the configuration of the robot whose deformation varies with time, the DOF vector is a function of time, i.e. q ≡ q(t). At time t = 0, the robot is undeformed and the DOF vector is q(0) ≡q; in the following,( ) represents evaluation of a quantity in the undeformed configuration. The strain calculations coming next all rely on the DOF vector q. Based on this kinematic representation, we will sequentially discuss the formulation of elastic energies and forces, external forces, and simulation loops.

Elastic energies and forces
The total elastic energy of a flagellated robot structure is the linear sum of stretching, bending, and twisting energies: where E s , E b , and E t are the stretching, bending, and twisting energies, respectively. Elastic energies of a structure are associated with the corresponding macroscopic strains, axial stretch, curvature, and twist [35]. The stretching energy associated with each edge is related to the axial stretch of the edge. Axial stretch is the change in length of an edge, normalized by the undeformed length. The axial stretch, k , of edge e k is where ē k is the undeformed edge length. Given this axial stretch, the stretching energy along edge e k is where E is the Young's modulus, A = πr 2 0 is the crosssectional area of the flagella, and r 0 is the flagellar radius.
The total stretching energy of the robot is the sum of individual stretching energies, i.e.
For edges on the rigid head and disk, the stretching stiffness EA is set to be sufficiently large to ensure that deformation is negligible. Bending energy is related to curvature, a node-based quantity that is related to the turning angle, ψ k , displayed in Fig. 3(c). Note that curvature is applied to all nodes except terminal nodes, and the curvature binomial for node x k is (κb) k defined as the following vector where (κb) k = 2 tan ψ k 2 . The scalar curvatures along the first and second material directors, calculated using the curvature binomial, are The bending energy is then calculated according to the equation indicates summation over all the discrete bending energies, ∆l k = 1 2 ē k−1 + ē k is the Voronoi length for the k-th node,κ (1) k andκ (2) k are the material curvatures in the undeformed configuration, and EI = π 4 Er 4 0 is the bending stiffness of the rod. For the rigid robotic head and 3D-printed plate, the value of bending stiffness EI is assumed to be so large that the curvatures at the rigid nodes remain nearly constant throughout the simulation.
Finally, the twisting energy is related to the relative rotation of the material frames between two adjacent edges, i.e. twist. The twist at the k-th node is where ∆m k,ref is the reference twist, which is the twist of the reference frame as it moves from the (k − 1)-th edge to the k-th edge. The method by which we calculate this reference twist is detailed at the end of Section 4.2 in our previous work [22]. The twisting energy is then calculated according to the equation whereτ k is the undeformed twist along the centerline, G is the shear modulus, and GJ = π 2 Gr 2 0 is the twisting stiffness. For the rigid components, the twisting stiffness is sufficiently large. Additionally, we assume that the material of flagella is nearly incompressible (i.e. Poisson's ratio ν = 0.5), so G = E/3.
In summary, at each DOF q k , the elastic forces (affiliated with nodal position) and elastic moments (affiliated with the twist angles) are where k = 0, 1, ..., ndof − 1.
In a single elastic rod, each internal node is associated with a discrete bending and twisting energy. However, in this paper, the flagellated robot is represented as a network of rods and a "T"-shape joint node. This node is associated with multiple discrete bending and twisting energies(denoted as d.B.E in Fig. 3(b)). In order to simulate the dynamics of the robot, EOM for each DOF is required, which includes not only elastic forces but also external forces. Consequently, we illustrate how we model the external forces from the GM onto the robot in the next section.

Fully-implicit simulation
In order to simulate the locomotion of the robot, time is discretized into small time steps of size ∆t. At each time step t j , the DOF vector q is updated. Using the following equations of motion (EOM), the k-th DOF marches from t = t j to t = t j+1 = t j + ∆t: where q k (t j ) andq k (t j ) are the known DOF and velocities at the previous time step, respectively, E E is the elastic energy evaluated at t j+1 , F ext k is the external force (or torque for twist angles) on the k-th DOF, and m k is the lumped mass at the DOF. The external viscous force exerted by the GM onto the flagella are detailed as Eqs. 2-4 in Section III A 1.The external viscous force and torque onto the head are calculated through Eqs. 23 and 24. Only the subscript k needs to be added to denote the force/torque at k-th node. Because the dynamics of this system is dominated by viscous forces with negligible influence of inertia, the results presented are not mass dependent as long as low Reynolds number is maintained. Eq. 40 represents the collection of ndof equations that has to be solved to get the new DOF q k (t j+1 ). Essentially, this equation is a statement of "mass times acceleration = elastic force + external force" at the k-th DOF. Once q k (t j+1 ) is updated, the velocity at time t j+1 is determined asq k (t j+1 ) = (q k (t j+1 ) − q k (t j )) /∆t. The EOM are solved using the Newton-Raphson method, i.e.
where f is a vector of size ndof, the k-th component of this vector can be computed from Eq. 40, and J is a square matrix representing the Jacobian for Eq. 40. The (k, i)-th element in the Jacobian matrix is where δ ki represents Kronecker delta, the terms are gradient of inertia, elastic forces, and external forces, respectively, in the order shown in Fig. 40. Well-documented evaluation of the gradient of the elastic energy ( ∂EE ∂q k ) as well as its Hessian ( ∂ 2 EE ∂q k ∂qi ) can be found in Jawed et al. [34] and Bergou et al. [36].

Remarks on algorithm
Next, we summarize the novelty of our algorithm for simulating a robot with multiple flagella described above. Note that solving Eq. 41 is the most computationally expensive part of the entire simulation procedure. It is crucial to notice the sparsity of the Jacobian matrix, J, and exploit its sparsity during the solution process [37], which helps reduce the computation cost. If the structure to be simulated is a single elastic rod (unlike a network of rods in this paper), the Jacobian is banded, and the time complexity of this algorithm is O(N ) [36]. However, the Jacobian in this paper is not banded due to the presence of the joint node x 2 in Figs. 3(a)(c). Referring to Fig. 3(b), the entire structure is a combination of stretching springs (e.g. one stretching spring is between x k and x k+1 ) and bending-twisting springs (e.g. one bending-twisting spring is between x k−1 , x k , and x k+1 ). The stretching energy of each spring (Eq. 32) only depends on six DOFs (nodal coordinates of two nodes). For the stretching spring on edge e k , these DOFs are x k and x k+1 . The gradient vector ∂ ∂q 1 2 EA k 2 ē k has only six non-zero terms and the Hessian matrix ∂ 2 ∂q∂q 1 2 EA k 2 ē k has only 6 × 6 non-zero terms. As for the bending and the twisting energies of each spring (Eqs. 36 -38), they are only dependent on eleven DOFs, i.e. x k−1 , θ k−1 , x k , θ k , and x k+1 in case of the spring at x k in Fig. 3(b). As a result, the gradient vector and Hessian matrix corresponding to these two energies contain only eleven and 11 × 11 non-zero terms, respectively. The complete expressions for the gradient and Hessian terms are available in works [34,38,39]; the coding implementation can also be found in open-source repositories [39][40][41].
Last but not least, a crucial contribution of this study is the observation that the actuation of the robot, e.g., the rotation of the motor, can be readily accounted for in the framework above by updating the undeformed configurations with time. Normally, the undeformed configuration of a structure is fixed and assumed to remain constant throughout the simulation. The strains in undeformed configuration (e.g.κ (1) k ,κ (2) k ,τ k in Eqs. 35a, 35b, and 38) are used in calculation of elastic energies, their gradient (i.e. elastic forces), and Hessian. However, in this class of robots, the rotation of the motor causes the undeformed twist at the head node (x 1 ) to vary with time. At each time step t j , the robot is actuated by updating the undeformed twist of the first node according to the rotational speed of the motor, ω T : This actuation causes the rotations of the head (ω h ) and the flagella (ω t ) along opposite directions such that |ω T | = |ω h | + |ω t |. The total rotational speed, ω T , is a control parameter in this study while the other one is the number of flagella, n.

C. Physical parameters
The following session will illustrate how we vary the robot design and experimental setup to realize designs 1 and 2. Remember that design 1 corresponds to the phenomenon that the propulsive speed of the robot decreases as the number of flagella increases, whereas design 2 reverses this phenomenon. Because the robot designs utilized in these two designs are different, we will provide the geometric and material parameters of the robots during experiments as follows. Hereafter, if we do not specify which design the parameters apply to, they apply to both. The number of flagella is n = 2 or 3 and each flagellum has a length of L = 0.111 m (design 1) and L = 0.089 m (design 2), a radius of r 0 = 3.2 mm, a density of 1000kg/m 3 (which is used to compute m i in Eq. 40), a Young' modulus E = 1.2 × 10 6 N/m 2 , and Poisson' ratio ν = 0.5 (incompressible material). Referring to Fig. 3, the radius of the robot head is R h = 0.02 m (design 1) and R h = 0.015 m (design 2), and the diameter of 3D-printed circular plate is L p = 0.04 m (design 1) and L p = 0.03 m (design 2). Recall that parameters C 1 , C 2 and µ in Eqs. 25 and 27 are fitting parameters and will be fitted later in Section V C 1. The time step used in this paper is ∆t = 10 −2 s and the length of each edge on flagella (in undeformed state) is ē k = 4.11 mm. Convergence studies were performed to ensure that the size of temporal and spatial discretization (∆t, ē j ) has negligible impact on the simulation performance.

IV. EXPERIMENTAL DESIGN
As the flagellated locomotion is an intricate interplay between the elasticity of the flagella and hydrodynamic loading from the surrounding GM, our experiments considered both the robot design and selection of GM for propulsion. The robot is placed inside a transparent cylindrical tube filled with GM as shown in Fig. 4 (a). Before finalizing the cylindrical tube as the container, we initially placed the robot at the midline of the height in a rectangular tank of 50×50×40cm or a circular tube with a radius of 25mm, but the robot moved closer and closer to the surface of the GM because of a drag-induced lift [30] though the robot was manufactured as symmetrically as possible.

A. Robot structure
Given that Euler-Bernoulli beam theory in Section III A predicts the presence of designs 1 and 2, experiments were performed to realize the prediction. Similar to the compact and lightweight robot used to realize design 1 (see Fig. 2 in our prior work [22]), the photograph in Fig. 4(b) depicts our robot for design 2 realization, which is propelled by n soft elastic flagella. It consists of four distinct components: multiple straight elastic flagella, Fig. 4(c) one 3D-printed circular plate attached to the motor shaft that protrudes from the robot head to hold flagella, and a head embedded with Fig. 4((d)) eight steel bearing balls (G25 Chrome Steel-AISI 52100) with an approximate diameter of 7.74 mm and mass of 2.06 grams, Fig. 4((e)) two 3.7V 200mAh rechargeable 502025 LiPo batteries (from Du litter energy battery), and Fig. 4((f )) one DC geared motor (from uxcell) with a nominal voltage of 3V, a power output of 0.35W, and a stall current of 0.55A. Inspired by natural creatures such as scorpions, snakes, and sand lizards, which have heads that are not spherical but more triangular, the head of our robot is designed as a cuboid (see Fig. 2 in Du et al.) or nose cone (see Fig. 4(b)). According to Euler-Bernoulli beam theory-based analyses, illustrated in Section III A 2, a smaller ω h is advantageous for the realization of design 2. As a result, compared to the previous cuboid head design, the head in Fig. 4(b)) has not only a different shape but also higher inertia and a more frictional surface. Metal balls in Fig. 4((d) are added to increase the inertia. Both head shapes contribute to the fluidization of the GM in front of the robot. To fabricate the elastic flagella, we applied the molding and casting techniques developed by Lazarus et al. [42] and Miller et al. [43]. The silicone-based rubber (vinylpolysiloxane from Elite Zhermack) injected into a PVC tube (e.g. from VWR International) mold and the inner and outer diameters of the PVC tube are the same as in our previous work [22], 3.175 and 6.35 mm. The radius of elastic flagella can be varied by using PVC tubes with different inner and outer diameters, allowing our robot platform to be effortlessly scaled up or down. The relationship between n and robot speed v is investigated in detail and will be explained in Section V A.
The spinning soft flagella propel the robot at a rotational speed of ω t . Since the torque is system-balanced, the head is then actuated at a rotational speed of −ω h . The control parameter is the rotational speed of the flagella relative to the head, ω T = ω t + ω h . To vary ω T , we replaced the geared motors while maintaining the other components the same instead of adding an encoder, which would significantly increase the robot size. The size and weight of all the motors are almost the same, (15-17) × 12 × 10 mm (L×W×H) and 13-15g even if they supply different ω h . When necessary, we added electrical tapes around the motor to account for the small differences in size and weight for different motors. Moreover, the value of ω t is proportional to the battery voltage and will decrease when the voltage drops. Subsequently, to maintain a constant ω t , we fully charged the batteries before each experiment and recharged them after each approximately ten-minute experimental trial.

B. Granular medium
A granular medium has to be chosen to complement the untethered flagellated robot design. Too much internal friction between granules prevents the robot head from fluidizing the granules in front and propelling forward, while too little friction results in insufficient flagellar propulsion. Internal friction during locomotion is correlated with the phase transition of GM from solid to fluid, which is controlled by temperature and volume fraction (VF) [44]. The less effort required to alter the VF of GM, the more beneficial it is for studying the locomotion in GM. Apart from the reason mentioned above, we chose gel soil water crystal beads (from EBOOT) as the GM for their transparency. Robotic performance, e.g., the robot position and the rotational speeds of head and flagella, could be viewed outside the medium using a traditional digital camera (Nikon D3400). The beads had a diameter of 2.5 mm when dried and d f = 9.4 ± 0.4 mm when completely saturated with water. The size of beads is controllable by adjusting the time they are submerged in water and is reversible after dehydration. During experiments, the robot for design 1 ran in the GM with d f = 9.4 ± 0.4 mm whereas the one for design 2 ran in the GM with d f = 5.2 ± 0.4 mm (see Fig. 4(b)). The corresponding VF values were approximately 0.52 and 0.54. In the former case, v decreases as n decreases (design 1), while in the latter case, the converse (design 2) is true. Throughout all experiments in designs 1 or 2, the temperature variation of GM is kept to a maximum of 0.5 degrees, and the granule configuration, including size, density, and homogeneity, is maintained.

C. Locomotion in granular medium
We used a transparent acrylic cylindrical tube (from FixtureDisplays, Amazon) with an inner radius of 53 mm and an axial length of 1220 mm as the reservoir for the GM to conduct the locomotion experiments. The tube was filled with GM and placed horizontally, perpendicular to the direction of gravity. Before each experiment trial, we placed the robot at one end at the center of the cross-section of the tube. Since the force applied by the GM on the robot in the cylindrical tube is axially symmetric, the confinement effect of the relatively close bounding wall is canceled out. As the robot was positioned at the center and thus surrounded by compact granules, the drag-induced lift [30] was suppressed. Hence, the robot actuated by the rotating flagella would move in GM along a roughly straight line. The video camera captured this movement at a frame rate of 29.98 fps. In addition, both the head and flagella were marked with markers of a different color than the corresponding robot components. For instance, as shown in Fig. 4(b), a black marker was attached to the yellow-colored head, and a red marker was attached to one of the dark green colored elastic flagella. Together with the GM's transparency, this operation ensured the accuracy of counting the rotational speeds of the robot head (ω h ) and flagella (ω t ) and locating the position of the robot, s in recorded experimental videos.

V. RESULTS AND DISCUSSION
Recall that there are two control parameters in our study, the total rotational speed of the motor embedded in the head, ω T , and the number of flagella, n. As shown in Fig. 1, when the motor is powered on, the actuation from the motor, ω T , is split into a constant rotational speed of the head, ω h , and flagella, ω t , and the rela-tionship ω T = ω h + ω t is always satisfied (ω T , ω h , and ω t are all non-negative values). Because the system is torque-balanced, the head and flagella rotate in opposite directions. The net propulsive force is the residual of the propulsion generated by flagellar deformation (with the deflection w) and the drag on the head and flagella from the GM, propelling the entire robot forward (along with −y direction) at speed v as illustrated in Fig. 2. In this section, we display the result comparison among the experiments and beam theory-based and DER-based numerical simulations. Since there are two designs, the sections below cover two sets of fitting parameters, experimental data, and numerical results.

A. Speed of the robot
The speed of the robot is the rate of change of its position with respect to the corresponding time interval, i.e. v = ∆s ∆t . Except in some cases where sometimes the position of the robot remains unchanged (and "jamming" happens) when the number of flagella is large (e.g. n = 4 and n = 5), the position of the robot increases proportionally with time. Hence, v is a constant throughout each experiment trial. Since we apply RFT to model the hydrodynamic forces onto the robot from the GM, periodic "jamming" phenomenon cannot be captured by our simulator. We focus only on robots with n = 2 and n = 3 flagella that maintain a constant v over time. However, "jamming" (i.e. the viscosity of the GM in front of the head increases) can be effortlessly modeled as a function of the robot configuration and integrated into the DER-based simulator. We exhibit experimental results in Figs. 5(a)(b) (design 1) and (c)(d)(e)(design 2). To obtain each data point ((ω T ,v) in Fig. 5(a)(e) and (ω T ,ω h ) in Fig. 5(b)(d) wherev = vη p L 4 /(EI) andω h = ω h η p L 3 /(EI) ), we randomly selected three separate one-minute sequence from every ten-minute experimental trial, repeated the operation for three different experimental trials, processed the data for every oneminute sequence, and calculated the average and variance (of nine values). Notice that there are no data points betweenω T = 0 andω T ≈ 200 in Figs. 5(a)(b) andω T = 0 andω T ≈ 18 in Figs. 5(c)(d). This is because we observed from experiments that if ω T is below a threshold (ω s T ), the robot would stay stationary (v = 0) but started to move continuously otherwise. The threshold changed depending on the robot design, such as the shape, surface friction and inertia of the robot head and flagellar length. For example, ω s T ≈ 50 rpm (ω T ≈ 200) for the robot used in design 1 [22] while ω s T ≈ 33 rpm (ω T ≈ 18) for the robot in design 2.
Note that the robots for designs 1 and 2 have distinct characteristics, such as the length of flagella L, and all data in Fig. 5 normalized and hence dimensionless. Normalization involves design parameters, L for example. As a result, though the ranges ofω T in Fig. 5(a)(b) and (c)(d) are seemingly different, they really overlap quite a bit (ω T = 0 − 250 rpm in Fig. 5(a)(b) and ω T = 0 − 100 rpm (Fig. 5c)(d)). To keep the same number of data points, and make ω T low enough for Euler-Bernoulli beam theory to apply, we did not include any additional data in Fig. 5(c)(d).

B. Parameter fitting for simulations
Next, we will show the DER-based and beam-based numerical simulation results for locomotion modeling. Recall that C 1 (in Eqs. 25 and 23) and C 2 (in Eqs. 27 and 24) are two fitting parameters to account for the shape and surface roughness. Additionally, µ, the constant used to quantify the body-granule friction coefficient, is the third fitting parameter. As detailed next, C 1 , C 2 , and µ in both designs 1 and 2 were generated from experimental data of 2-flagellar locomotion, and the same values were used in simulations to predict 3-flagellar locomotion. The predicted outcomes were compared to the corresponding experimental data and validated the fitting process.
The normalized total rotational speed of motor,ω T , compares the period of angular rotation to the elastoviscous relaxation time,T = η p L 4 /(EI), where EI is the bending modulus of the flagella and η p is the drag coefficient along the perpendicular direction [46], i.e. ω T = ω T η p L 4 /(EI). Similarly, the rotational speed of the head is normalized asω h = ω h η p L 4 /(EI) and the velocity of the robot is normalized asv = vη p L 3 /(EI). In Figs. 5(a)(b) (for design 1), we plot the normalized speed of the robot,v, versus the normalized total rotational speed of the motor,ω T and normalized rotational speed of the headω h versusω T , respectively. All the experimental data, DER-based and beam-based simulation results for n = 2 and n = 3 are shown in the figures. The experimental data (v vs.ω T andω h vs.ω T ) for n = 2 were adopted to determine the best fit values of the fitting parameters: C 1 = 2.420, C 2 = 0.039 and µ = 6.828 (for design 1). These parameter values were then applied into the numerical simulators to simulate the locomotion performance of n = 3. In design 2, the same technique was applied in Figs. 5(c)(d) to find the fitting parameters that minimize the fitting error between experimental and DER-based simulation results : C 1 = 28.750, C 2 = 0.938 and µ = 2.125 (for design 2). Figs 5(a)(c) demonstrate a high degree of agreement between experiments and DERbased simulations. The robot designs for designs 1 and 2 are distinct, and the relevant physical parameters were detailed in Section III C.
The following modeling assumptions can partially cause the slight disparity between the experimental and simulation results: • In this fluid (i.e., the granular medium) structure (i.e., the robot) interaction modeling of robot locomotion, we assume that RFT can characterize the drag from the GM (the fluid).
• The structure model (DER) assumes that the flag- ella are infinitesimally thin elastic rods.
• It is assumed that the drag force acting on the head is considered to be linearly proportional to the velocity, and torque on the head is linearly proportional to its angular speed.
In addition, inevitably, there are experimental errors, structural defects introduced during fabrication, for instance. Nevertheless, the reasonable consistency between experiments and simulations suggests the validity of RFT in this context.

C. Speed vs. robot geometry
To quantify the effect of geometrical design, e.g. the number n and length L of flagella, we perform a parameter sweep along the angular velocity to systematically study the mechanical response of robots in both designs.

Speed vs. number of flagella
After performing the first round of experiments and data processing, we observed a counter intuitive phenomenon from Fig. 5(a) that the robot with two flagella (n = 2) moved faster than the one with three flagella (n = 3), at the same value ofω T . Moreover, the relationship betweenv versusω T is nonlinear. We then used the optimal set of fitting parameters described in Section V B to perform the beam-based framework calculation; the results ((v vs.ω T andω h vs.ω T )) are plotted in Figs. 5(a)(b). We denote points A and B in Fig. 5(a) as the intersection of beam-based and DERbased simulation results for n = 2 and n = 3 in design 1, respectively. As can be observed, the beam-based simulation results are close to DER-based simulation results whenω T is less than the value at point A (n = 2) and point B (n = 3), but deviate in other cases. All these observations indicate the significant flagellar deformation and the tight coupling between the head and flagella.
Nonetheless, the beam-based mechanics analysis framework in Section III A 1 tells us that if there is no head or the head is fixed along the x-axis (in Fig. 2), F p will increase proportionally to the increase in n. The propulsive speed of the robot, as illustrated in Fig. 2, is the result of a complex battle between the projection of v t along −x (propulsion) and the sum of F p in Eq. 23 and the projection of v p along −x (friction force). Section III A 2, on the other hand, elucidates that the effect induced by the head is one of the reasons for the difference between designs 1 and 2. Notice thatω h occupies more than 80% ofω T in Fig. 5(b), indicating the majority of the increment inω T is spent growingω h . As a result, to achieve design 2, which is predicted by the beam-based framework, we have to slow down the rotational speed of the head. As shown in Fig. 5(d), the ratio ofω h toω T in design 2 drops to around 55%. Section IV A details how this was accomplished by updating the robot design in design 2 relative to the one utilized in design 1. We performed experiments with the updated robot design and the results are shown in Figs. 5(c)(d).
Note that all of the results that pertain to beam theory in Fig. 5 are NLB w/o head, as illustrated in Section III A 2. As mentioned in the last session, when n equals two, the parameter fitting technique was leveraged to match the experiments and simulations. Point C in Fig. 5(c) stands for the conjunction of DER-based simulation results of n = 2 and n = 3. Experimental data point out that a two-flagellar robot moves slower than a three-flagellar robot while all other parameters remain constant. There is no evident difference inv between n = 2 and n = 3 simulation results for values of ω T smaller than the value at point C, which was captured by experiments. However, beam-based prediction does not show design 2 in Fig. 5(c) but it demonstrates design 2 when we increase the range ofω T as shown in Fig. 5(e). Furthermore, our goal of reducing the proportion ofω h taken fromω T is verified in Fig. 5(d). In Fig. 5(e), we display the beam-based framework predictions with the fitting parameters C 1 = 28.750, C 2 = 0.938 and µ = 2.125 plugged in. Here, D represents the watershed point, beyond which the robot with n = 3 moves faster than the robot with n = 2 and vice versa. The beam-based prediction follows the same trend with the DER-based simulation in Fig. 5(c) even though the watershed points are drastically different. Also, the magnitudes ofv andω h predicted by the beam-based framework are a lot larger than those obtained from experiments and DER-based simulation. This further supports our assertion that robot locomotion involves large flagellar deformation and close interaction between the head and flagella. In summary, the simple beam-based analytical framework qualitatively captures the relationship between the speed of the robot and the number of flagella, though not quantitatively.
One more thing that needs our attention is the nonlinear increment of slope in Fig. 5(a). Both slopes for n = 2 and n = 3 initially increase in magnitude until they reach a point after which they begin to decrease in magnitude. We refer to the point for n = 2 as point E and the one for n = 3 as point F . The points E and F are conspicuous because, at these points, the unit magnitude increase inω T turns into the steepest increase in v, indicating the maximum efficiency of the robot locomotion is reached. The same phenomenon happens in design 2 (see Fig. 5(c)). Last but not least, point E is near point A and F is in the neighborhood of point B. This could mean that an increase in large flagellar deformation implies a decline in locomotion efficiency.

Speed vs. flagellar length
As previously stated, the data in Fig. 5 has been normalized. To help readers visualize the physical scenarios, we plot the simulation results from DER and NLB in Fig. 6, which depicts the relationship between the speed of the robot v in terms of each flagellar length, L, and the total rotational speed of the motor, ω T . As illustrated in Fig. 6(a)(design 1), DER demonstrates that when the number of flagella, i.e. n is fixed, the distinction between DER and nonlinear beam theory (NLB ) becomes more discernible as flagella become longer. When ω T is small, DER results indicate that v grows as the value of L increases, but this trend breaks down when ω T is big. In comparison, NLB behaves more consistently and linearly, i.e. when n is constant, v increases with the growth in L. These observations make sense in light of the fact that we previously acknowledged that DER is capable of capturing the nonlinearity of flagellar deformation, whereas NLB is not. The DER results displayed in Fig. 6(b)(design 2) show that as L increases, the difference between v of n = 3 and n = 2 does not always increase, but in the range of ω T presented, v of n = 3 is greater than v of n = 2. In comparison to Fig. 5(c), NLB predicts the appearance of design 2 after a threshold is reached. However, this threshold becomes smaller when L is larger. In reality, we tested designs 1 and 2 on robots with n = 4 and n = 5 flagella. Nonetheless, "jamming" happened from time to time randomly. As a result, only simulation data are presented here.

D. Deflection of beam end
In the last session, we state that the simulation results from the beam-based framework (NLB w/o head design in Section III A 2) are qualitatively compatible with both experimental and DER-based simulation results. We also emphasize the importance of substantial flagellar deformation in relation to the performance of robot locomotion. Consequently, in Fig. 7, we compare beam-based deflection in three regimes, LB, NLB, and NLB w/o head as described in Section III A 2. The simulation parameters are identical to those for design 1 in Section III C and to the fitting parameters for design 1 in Section V B. As seen in Eq. 18, the effect of the head The plot of robot movement speedv versus total rotational speed of the robotωT ; (a) corresponds to design 1, while (b) corresponds to design 2. In the abbreviated legend, "D" stands for DER-related results and "B" represents results from the nonlinear beam theory-based framework; the number following "D"/"B" represents the length of each flagellum in millimeters(mm), "n2" denotes two flagella and "n3" denotes three flagella. For example, "D130n2" stands for the results of a DER-based simulation of a robot with two 130mm flagella.
in the flagellar deflection is ignored in regime LB, so the rotational speed of the flagella, ω t is chosen as the independent variable. To see more noticeable difference in flagellar deflection among three regimes, we picked the scenario with larger flagellar deformation, design 1, and the corresponding fitting parameters to compute ω h in NLB and NLB w/o head regimes. As seen in Fig. 7, when ω t is relatively small (ω t 5 rpm andω T 20), the flagellar deformation drops to the LB regime introduced in Section III A 2 and the difference in deflection of the end of the beam (flagella) among the three regimes is subtle ( 5%). Moreover, as can be observed in Fig. 7, the deflection of the end of the beam in regime NLB w/o head is almost always equal to that in NLB throughout the range ofω T . This corroborates our statement that the head velocity-induced effect could be ignored in Section III A 2.

VI. CONCLUSION
In this study, we designed a low-cost experimental setup, an untethered robot actuated by multiple soft flagella moves in granules. Meanwhile, we developed two numerical simulators, one based on Euler-Bernoulli beams and the other on discrete differential geometry (DDG) framework, to simulate the performance of the articulated locomotion in granular media (GM). Both numerical tools use resistive force theory (RFT) to model the drag force exerted by the GM on the robot and Stokes' law to model the external force/moment applied to the robot head by the GM. Initially, experiments unveiled a counterintuitive phenomenon: the robot's speed decreases as the number of flagella increases (design 1). However, the beam-based simulator predicts the existence of the converse case (design 2) on the condition that the rotational speed of the head is suppressed to be zero. The robot design was then modified to slow down the rotational speed of the robot head. The DDG-based simulator models the robot into a composition of Kirchhoff elastic rod and discretizes it into a series of mass-spring systems. As a result, the elastic energy of the robot structure is the linear sum of the discrete elastic energies associated with each spring, which include the stretching and coupled bending and twist-ing energies. At each time step, the equations of motion formulated can be summarized as follows: at each degree of freedom (DOF), the sum of elastic forces (i.e., the negative gradient of elastic energies) and external forces equals the lumped mass multiplied by the acceleration of that DOF. Especially, the actuation of the robot (i.e., the rotational speed of the motor) is modeled by a timevarying natural strain (twist) at the node standing for the head. This method enables us to simulate the robot locomotion fully implicitly.
Both beam-based and DDG-based frameworks capture designs 1 and 2 successfully. Although the beambased simulator can quantify only the trend in both designs while the DDG-based method accurately reproduces experiments, it is still exhilarating because the simple beam theory-based analysis depicts the complicated locomotion system. Simultaneously, the discrepancy between beam-based prediction and experiments indicates the large flagellar deformation and nontrivial coupling between the head and flagella. Additionally, both simulators, validated experimentally, shed light on the highly nonlinear functional relationship between the locomotion performance, such as speed and efficiency, and its physical parameters, such as the number of flagella. The nonlinear dependency of the speed of the robot concerning the rotational speed of the motor necessitates the development of a design tool for optimal control of this class of robots. Last but not least, thanks to the simplicity of the beam-based analysis framework and the computational efficiency of the DDG-based simulator, they can be exploited to perform parametric studies and identify the optimal design and control this class of articulated robots as long as they locomote through the GM. We chose a smooth and soft GM with a large particle diameter. Scaling power requirements for locomotions moving in harsher GM-like sand will be an interesting challenge.