Gianmarco Mengaldo

Federico Renda
email: federico.renda@ku.ac.ae

Steven L Brunton
email: sbrunton@uw.edu

Moritz Bächer
email: moritz.baecher@disneyresearch.com

Marcello Calisti
email: mcalisti@lincoln.ac.uk

Christian Duriez
email: christian.duriez@inria.fr

Gregory S Chirikjian

Cecilia Laschi

Renda Federico

A concise guide to modelling the physics of embodied intelligence in soft robotics

Keywords: Soft robotics, Embodied intelligence, Mathematical modeling, Computational modeling, Computational physics

établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Soft robotics is largely motivated by the functional role of soft tissues in living organisms [START_REF] Kim | Soft robotics: a bioinspired evolution in robotics[END_REF]. Life has had millions of years to adapt to its surrounding environment and co-evolve nervous and muscle-skeletal systems to achieve task-efficient performance, synergistically. We observe that living beings are soft and compliant, and we argue that this is instrumental to their embodied intelligence [START_REF] Pfeifer | Self-organization, embodiment, and biologically inspired robotics[END_REF]. According to this modern view of intelligence, the physical body play a much larger role in shaping intelligence, since a part of sensory-motor behaviour emerges from its interaction with the surrounding environment, with minimal or no involvement of the nervous system. Soft biological systems use their complex internal body structure to efficiently leverage physical interactions with the external environment and achieve the desired actions. Indeed, external interaction forces, instead of being treated as disturbances needing compensation, are used for the intended movements [START_REF] Blickhan | Intelligence by mechanics[END_REF]. As an example, in locomotion, gravity is exploited for stepping forward, and adaptation to uneven terrains is provided by compliant elements within the leg joints, with limited need for active inputs from the central nervous system. Similarly, octopuses, an iconic model for soft robotics, adopt highly effective unfolding arm reaching movements by leveraging the buoyancy and interaction forces from the water surrounding them.

Embodied intelligence is a well-known paradigm in robotics and in soft robotics. In a typical robot sensory-motor behaviour scheme, we can consider a sensory system perceiving the environment, controllers processing the incoming information and planning a motor action, and then a mechanical system that executes motor actions in the physical environment. Embodied intelligence can be seen as the mechanical feedback from the environment, directly onto the mechanical system of the robot physical body, with no involvement of controllers or processing (see Fig. 1, from [START_REF] Pfeifer | Self-organization, embodiment, and biologically inspired robotics[END_REF]). This view gives a clear perception of how powerful embodied intelligence can be in simplifying robot sensory-motor behaviour and increasing overall robot efficiency and effectiveness. All this works if we assume that a soft, compliant body receives, and does not reject, the mechanical feedback from the environment. Fig. 1: In a typical sensory-motor scheme, either for a robot or a biological system, three main components are usually considered: a sensory system (receptors in biology) receiving inputs from the external environment, some form of controller (nervous system in animals), and a mechanical system (musculoskeletal system in biology) acting on the environment for accomplishing a task. Embodied intelligence can be seen as the mechanical feedback received by the physical body from the environment. It involves the body dynamics and morphology with no involvement of the controller. It allows closing a very short control loop, by-passing most of the computational processes. (Reprinted with permission from [START_REF] Pfeifer | Self-organization, embodiment, and biologically inspired robotics[END_REF] and adapted).

How to systematically design embodied intelligence into soft robots is one of the current challenges in the field of soft robotics. We argue that this goal can be achieved by means of a mathematical description of the physical interactions characterizing embodied intelligence. Related computational modeling would enable simulations to be used for both design and control purposes.

The task of constructing this modeling framework is not an easy one. We focus on the soft parts that are necessary to implement embodied intelligence and that go beyond the rigid-body modeling techniques used in robotics. Take as an example a soft arm immersed in a fluid and interacting with solids (see, e.g., the octopus in Fig. 2). To describe and model the physics underlying embodied intelligence, it is necessary to consider:

• the continuous deformations of the arm deriving from muscle activations,

• the coupled interaction of the arm with the fluid, when e.g., the arm reaches out for target objects and moves the fluid , and the fluid contributes to its deformation,

• the interactions of the soft arm with solids (rigid or deformable), like the seafloor when walking, or external objects when grasping.

The description of these three points present a number of significant challenges.

First, the problem is inherently multiphysics, due to the physically heterogeneous nature of interactions underlying embodied intelligence. In fact, one needs to take into account the physics of the soft body and the dynamics of muscle activation, the flow physics and its coupling with the soft body, and the physics describing contact and adhesion between two solids. This poses the critical challenge (and opportunity) of interdisciplinary work across multiple communities, ranging from robotics, fluid dynamics, structural mechanics, tribology, and contact mechanics. These fields typically have their own customs and terminology, and cross-fertilization can prove difficult, yet will be beneficial to bring advances in modeling embodied intelligence.

Second, the problem is multiscale, as the number of scales one needs to describe may range from millimeters to meters. For instance, in our example the arm reaching movement may undergo an overall displacement of several centimeters. Yet, the description of the deformation of the soft body, and its interaction with the flow and solid, may require a much finer description to accurately capture its behaviour. To this end, fast computational methods to solve multiscale problems are required. Integration with the scientific computing and applied mathematics communities may therefore prove important.

Multiphysics and multiscale problems are notoriously challenging but have successfully been conquered in some fields, including the aerospace industry, biomedical engineering, and material science, to cite a few. In these fields, it may be possible to develop a digital twin, which is a virtual representation of the system comprised of a hierarchy of models for different components and levels of fidelity, ultimately pinned to data [START_REF] Niederer | Scaling digital twins from the artisanal to the industrial[END_REF][START_REF] Brunton | Data-driven aerospace engineering: reframing the industry with machine learning[END_REF]. In the context of soft robotics, these problems have yet to be addressed, although promising advances have been made. Some of the factors that make modelling particularly difficult here are: i) how the actuation forces are designed and mathematically modelled to achieve a desired action, and how the surrounding environment influences them; ii) large deformations of nonlinear materials and their coupled interactions with the surrounding medium, which lead to high-dimensional models that are often prohibitively expensive to simulate and use inside control loops; and iii) partially known interfacial physics and unmodeled dynamics -e.g., nonlinear friction mechanisms for solid-solid interaction, and turbulent drag for fluid-structure interaction. In addition, the range of morphologies (e.g., arms, fingers, legs, fins), materials (e.g., hyperelastic, heterogeneous, functional), abilities (e.g., reaching, grasping, walking, morphing, growing, swimming, jumping, crawling, digging) and intended applications (e.g., healthcare, manufacturing, underwater sensing and manipulation, scientific exploration, entertainment) is extremely diverse, exacerbating the complexity of the modelling task. The latter aspect may lead to staggered and application-specific modelling strategies, that may not be beneficial to the soft robotics field.

In this paper, to support our argument of distilling embodied intelligence into physical interactions, we provide a concise guide to the latter and to their underlying mathematical models. In particular, we present and discuss the most prominent models that describe 1) the interactions of soft robot components (including actuators) that lead to efficient movements and deformation, and 2) the robot interactions with the surrounding environment. We refer to the former as internal interactions, and the latter as external interactions, where the environment can be constituted by a fluid or a solid. 1 To introduce the models underlying the physical interactions found in embodied intelligence, we use the octopus depicted in Fig. 2 as a proxy for a general biological or soft robotics system. The octopus represents an excellent model for soft-body parts, while not excluding rigid components in the overall treatment of the problem. Analogously, we consider interaction with a general fluid, which holds for liquids as well as gases. The description of the modeling techniques has the aim of providing the basic tools for their uptake. To this end, we give a comprehensive picture, including both internal and external interactions, both with solids and fluids. We also showcase their use and discuss the challenges that can promote beneficial interdisciplinary interactions. Fig. 2: An illustration of our scope in modeling the physics of embodied intelligence. We use an octopus as a proxy for a general soft body, in its environment, where we highlight the three key modeling areas that contribute to a substantial insight of embodied intelligence: internal interactions, and external interactions models, either with a fluid medium or a solid support. Internal interactions refer to the deformations of the soft body induced by actuators, or muscles. External interactions refer to the action of the external environment (fluid or solid) on the soft body. We also highlight how the accurate description of internal and external interactions can lead to low-dimensional models that capture the soft body behaviour with a reduced number of state variables.

A Mathematical Framework for Modeling of Soft Robots 150

The mathematical framework adopted belongs to the realm of mechanics that is at the foundations of physical interactions in both solids and fluids. To describe the mechanics of physical systems there are two views, Lagrangian and Eulerian. In the former, one tracks the trajectories of particles, while in the latter one observes the particle velocities at a fixed point in space. The Lagrangian setting is commonly used in classical solid mechanics, while the Eulerian setting is commonly used in fluid mechanics. These two views allow us to introduce an abstract mathematical formalism for the soft body subject to internal and external interactions that consists of a set of differential equations:

Dq sb = N sb + C int + C ext , in Ω sb . (1)
Equation ( 1), along with suitable initial and boundary conditions, can describe both models arising in continuum mechanics and multibody dynamics in the soft body domain Ω sb . The former typically yields a set of partial differential equations, while the latter a set of ordinary differential equations.

Depending on the approach used, the variables q sb , that describe the soft body (sb), can have a different meaning. For instance, in a Lagrangian view, q sb represents the position and momentum of material particles, while in an Eulerian view, q sb is the observed velocity at each given point. D is a differential operator that may be a partial (∂) or a total (d) derivative (including a material derivative D), of first or second order.

N sb is a nonlinear term that describes the soft body (sb) mechanics, which may depend on q sb , its partial derivatives with respect to spatial coordinates

x, time t, and by a set of tunable constants such as viscosity, stiffness, and other actuation and material parameters.

C int is a coupling term that accounts for internal (int) interactions (e.g., actuation forces for tendons or pressure of pneumatic chambers). C ext is a coupling term that accounts for external (ext) interactions (e.g., contact with external solids or interactions with the surrounding medium).

While equation (1) may not be the standard way of introducing models in the realm of soft robotics, it allows a sufficiently general mathematical framework that can be used as a basis to construct a unified formalism for the multiphysics of internal and external interactions. We divided the righthand side into a nonlinear term N sb , and two coupling terms C int and C ext , to emphasize the mechanics of the soft body (through N sb ) and its coupling not only with internal (through C int ) but also with external systems (through

C ext ).
In particular, C int , and C ext can constitute forcing terms, that results from lumping the interactions into simplified terms, or can be terms coupling the equation describing the soft body (1) to additional equations describing the physics of internal and/or external interactions (e.g., actuators, and the surrounding environment). These interactions could be for instance described by equality and inequality constraints, especially when the soft body is in contact with a medium. In this case, the constraints can be embedded into the equations by means of, e.g., Lagrange multipliers, and encapsulate the physics of interactions. For example, frictional contact between two surfaces is best modeled with inequality constraints which become active as soon as a soft body is in contact with either itself or the environment.

Starting from equation [START_REF] Kim | Soft robotics: a bioinspired evolution in robotics[END_REF], in the following we introduce the models and solution methods adopted for internal and external interactions. More specifically, the outline of this paper is as follows. In section 2, we detail the models for internal interactions, describing both the models for soft-body mechanics (i.e., N sb ), and the internal interactions with actuators (i.e., C int ). In section 3, we describe the models for external interactions (i.e., C ext ). In section 4, we outline how these models can be used in practice for soft robotics, such as by leveraging reduced order models. In section 5, we review the modeling challenges and opportunities described, and draw some conclusions and perspectives.

2 Internal interactions

Challenges

A soft-bodied system, such as the arm of the octopus depicted in Fig. 2, is composed of a set of sensors and actuators, distributed inside a soft tissue, that work together to achieve a given task. The various sensing and actuation mechanisms interact internally with one another, and they typically operate within a nonlinear and possibly heterogeneous material constituting the soft tissue. We focus here on the fundamental robotics task of modeling the relations between the actuator space and the soft-body deformations. A mathematical model describing internal interactions should therefore capture the robot deformation (i.e., N sb ) produced by the internal actuation forces, where the latter also needs to be accurately represented (i.e., C int ).

The physics of actuation corresponds to a subset of the total degrees of freedom of the system, associated to the physics of deformation of the structure. Once the subspace of actuation is defined, the key challenges specific to modeling internal interactions are:

1. to model the physics of the actuation in such a way that it can be identified on the real system, and 2. to couple the two physics: the one of the deformation of the structure, and the one of the actuation.

In terms of actuation, we take into account the following three main strategies Tendons (1D). We consider a punctual action of a cable or similar tendonlike actuator, that can be considered as one-dimensional (1D). If the tension is exerted by a tendon which pulls, it is relatively easy to describe the force field created on a soft robot structure, based on the geometrical path. When a series of tendons are placed on a robot and their input lengths need to be integrated, the problem becomes more complex: since the tendons are coupled, it is necessary to verify that they are tense to exert a force. The tension of one tendon on a structure could slack on other ones. Also the tension inside the tendon is a signed force: we can pull with a tendon, but not push. The algebraic equations, translating the motion of the tendon, lead to constraints of inequality and complementarity. Finally, more advanced models that could be used to replicate the mechanics of the tendons (extension, bending, internal friction, etc.) are neglected to focus on this tensile force.

Flexible fluidic actuators (2D). We consider chambers that deform when pressurized with a fluid, so that pressure acts on a surface, and is assumed to be constrained to two-dimensional (2D) deformations. In the same way, for these fluidic actuators, it is possible to directly take into account the pressure exerted in cavities. This pressure is then integrated on the surface of the cavity to calculate a force field applied to the soft robot structure. This approach is often used for pneumatic actuation, based on pressure regulators. On the other hand, in the case of liquid injection, it is often the variation of the cavity volume to constitute the input of the model. As with the tendon case, it is then necessary to add an algebraic equation in the model translating the volume variation. Moreover, the weight exerted by the liquid on the structure can create additional deformations. Finally, one can imagine, in the longer term, coupling the deformation models with dynamic models of the fluid.

Smart materials (3D). We consider a large class of materials responding to external stimuli that can be regarded as three-dimensional (3D) actuators. Some smart materials, such as electro-active or electro-ionic polymers, have their internal stress field dependent on an electric field. Other materials, such as shape memory materials, have their constitutive law dependent on temperature. Materials can also be integrated into soft robots that react to a magnetic field. In all these cases (and possibly in many others, the field of smart material being vast), the problem is the coupling of the deformation equations with other equations of physics, like electric or magnetic fields or temperature diffusion, often at scales that are much smaller than the scale of the robot.

The modeling of these actuators constitutes a challenge per se, and their integration with different soft-body models can be non-trivial. In particular, soft robots are highly under-actuated (the degrees of freedom associated with actuation are significantly fewer than the degrees of freedom of the soft body).

Therefore, it is frequently required to add algebraic constraints, often associated with Lagrange multipliers Λ and Θ, to drive the motion of the soft body within the actuator space. In this case, the coupling term with the actuators becomes C int = Λ T I e + Θ T I i , where I e and I i are equality and inequality constraints, respectively. The method of the Lagrange multipliers is a strategy to enforce equality constraints to a functional, in this case equation (1).

Models

Two modeling approaches of internal interactions exist in robotics. The first, called direct modeling approach, starts from the knowledge of the motion and/or forces in the actuator space and from other unactuated degrees of freedom to compute the robot shape. The second, known as inverse modeling approach, starts from a desired shape or position of the robot and calculates the actuator space in order to reach the desired shape/position.

In this section we will focus on the direct modeling approach. The actuator forces will modify the static equilibrium or the dynamics of the robot.

In the following, we introduce the main models that can be used for describing the soft body mechanics N sb and how they can be integrated with the actuation strategies reported above, through C int . We introduce continuum solid mechanics models (fully 3D or rod/shells) and finite-dimensional parametrization models and we keep their description at a level general enough to provide the reader with tools adaptable to most cases. We also introduce the possible use of data-driven methods for the same modelling problems.

Continuum three-dimensional solid mechanics models

The most general model that captures the full complexity of the soft body is given by three-dimensional continuum solid mechanics (see Fig. 3(a)). If we consider equation ( 1), taking an Eulerian view:

D becomes ∂/∂t, the partial derivative with respect to time t, q sb becomes v sb , representing the velocity field of the soft body, N sb becomes ∇ • σ sb , the divergence of the soft body stress tensor σ sb .

This equation describes the balance of linear momentum, that is usually complemented by the balance of angular momentum, through σ = σ T , and balance of mass, through ∂ρ/∂t + ρ(∇ • v sb ) = 0. One can also consider adding the balance of energy. The constitutive relations for describing the soft-body material usually rely on both elastic and hyperelastic models, or ad-hoc constitutive models. These ad-hoc constitutive relations can originate from e.g., phenomenological experiments of a specimen and provide an accurate simulation of the soft body deformation [START_REF] Xavier | Finite element modeling of soft fluidic actuators: Overview and recent developments[END_REF]. Note that hyperelastic models sometimes have many parameters that are difficult to identify in practice. Moreover, other phenomena such as viscosity, plasticity or anisotropy need also to be considered into the constitutive equations. Anisotropy has, for example, a direct influence on the kinematics of the soft robot [START_REF] Vanneste | Enabling the control of a new degree of freedom by using anisotropic material on a 6-dof parallel soft robot[END_REF]. Finally, these deformation models are very sensitive to the boundary conditions, thus to the modeling of the actuators but also to external interactions.

All the three actuation strategies, tendons, fluidic actuators, and smart materials, can be adopted in conjunction with continuum solid mechanics.

Their coupling with the soft body is obtained through the term C int . This can be constituted by Lagrange multipliers C int = Λ T I e + Θ T I i arising from a set of constraints imposed by the actuators on the soft body. They can also be imposed as a set of boundary conditions for equation [START_REF] Kim | Soft robotics: a bioinspired evolution in robotics[END_REF]. Solution strategies rely on numerical methods, e.g., finite and spectral element methods, whereby the partial differential equations are solved at a set of prescribed nodes (or mesh). This is typically extremely expensive computationally with a corresponding number of degrees of freedom in the order of O(N 3 ), where N is the number of nodes of the underlying mesh. Indeed, a convergent solution (i.e., a solution that provides the correct results) depends on the number and type of elements chosen as well as on their geometric quality, where complex geometrical features and nonlinearities might further influence the solution convergence. The requirements in terms of solution convergence might also differ depending on the way boundary conditions are applied -i.e., the convergence of one or another variable might be required. For instance, if we consider a fluid (liquid or gas) which inflates a cavity, and we load the cavity in volume (thanks to a constraint), then the convergence of the motion variables (position and velocity) is sufficient. However, if the structure is loaded in pressure, it must be in equilibrium with the stress variables of the material, hence, the convergence of stresses is needed. In addition, for the management of contacts with e.g., internal actuation systems, it is sometimes necessary to use a finer mesh on the surface of the objects in order not to alter their geometric definitions, thus undermining the accuracy of the simulation.

Successful implementations in the realm of soft robotics applications exist.

The software SOFA implements the finite element method to simulate continuum solid mechanics [START_REF] Duriez | Soft robot modeling, simulation and control in real-time[END_REF], and offers solutions, such as model reduction [START_REF] Goury | Fast, generic, and reliable control and simulation of soft robots using model order reduction[END_REF],

to find a compromise between accuracy and computation time. Fig. 4-(a)

shows an example of its use for modeling a soft pneumatic finger. ChainQueen implements a differentiable Lagrangian-Eulerian physical simulator based on the moving least square material point method for solid mechanics, along with actuation and contact with external objects [START_REF] Hu | Chainqueen: A real-time differentiable physical simulator for soft robotics[END_REF]. Evosoro uses a solid mechanics engine, Voxelyze [START_REF] Hiller | Dynamic simulation of soft multimaterial 3dprinted objects[END_REF], that allows the simulation of soft multimaterial robots [START_REF] Cheney | Unshackling evolution: evolving soft robots with multiple materials and a powerful generative encoding[END_REF]. These three options are open-source and currently under active development. Some additional commercial software to tackle the internal interaction modeling challenge exist, including ABAQUS [START_REF]ABAQUS software[END_REF], ANSYS [START_REF]ANSYS software[END_REF],

COMSOL [START_REF]COMSOL software[END_REF], and Altair [START_REF]Altair software[END_REF], among others. These provide platforms for solving full complexity solid mechanics problems, but they are not tailored to soft robotics.

This modeling approach, along with the computational strategies for solving it, allows the description of all topologies commonly required in soft robotics, e.g., rods, lattice of beams, shells, volumes. These topologies can be generated by means of Computer-Aided Design (CAD), and standard meshing practices. This however, may lead to an inexact geometrical representation of the soft body, that in turn can yield an inaccurate solution for the soft-body motion. This drawback could be addressed by using accurate mesh generation practices (e.g., high-order mesh generation [START_REF] Turner | Curvilinear mesh generation using a variational framework[END_REF][START_REF] Mengaldo | Industry-relevant implicit large-eddy simulation of a high-performance road car via spectral/hp element methods[END_REF]) and by increasing the number of mesh nodes, or with isogeometric analysis [START_REF] Cottrell | Isogeometric Analysis: Toward Integration of CAD and FEA[END_REF].

We finally note that the framework of 3D continuum mechanics allows for the description of robots with both soft and rigid components. To couple soft and rigid body components, soft-body degrees of freedom at the interface between the soft and rigid part can be expressed as rigid body position and orientation, or coupling constraints may be introduced [START_REF] Hoshyari | Vibrationminimizing motion retargeting for robotic characters[END_REF][START_REF] Adagolodjo | Coupling numerical deformable models in global and reduced coordinates for the simulation of the direct and the inverse kinematics of soft robots[END_REF].

Shell and rod models

A significant number of soft robots are characterized by an elongated structure with two dimensions much smaller compared to the third one. In this case, it is possible to greatly reduce the number of degrees of freedom required to describe the robot by adopting 1D rod models. Sometimes, only one dimension is negligible, and 2D shell models may be adopted. These fall under the realm of continuum mechanics although not in 3D.

A particularly useful continuum modeling approach to describe rods and shells in soft robotics is the Cosserat model. In this model, the 3D Boltzmann continuum with pointwise three (translational) degrees of freedom is replaced by a set of infinitesimal 6D micro-solids stacked along the dominant dimensions [START_REF] Antman | Nonlinear Problems of Elasticity[END_REF]. This intrinsic theory allows simulating the nonlinear finite rotation of rods and shells with better accuracy and efficiency compared to approaches deduced from a 1D or 2D reduction of the traditional 3D theory [START_REF] Meier | Geometrically exact finite element formulations for slender beams: Kirchhoff-love theory versus simo-reissner theory[END_REF]. For instance, 3D FEM corotational methods leads to a linearization of the kinematic equations that is found unsuitable for the analysis of arbitrarily large rigid-body displacements [START_REF] Shabana | Continuum-based geometry/analysis approach for flexible and soft robotic systems[END_REF].

The mathematical formalism underlying a Cosserat model relies on the Lie group of rigid body transformation SE(3), due to the assumptions made on the material microstructures. Following [START_REF] Boyer | Dynamics of continuum and soft robots: A strain parameterization based approach[END_REF], if we consider rods, a Cosserat rod is modeled by a continuous set of rigid cross sections stacked along a material line parameterized by a curvilinear coordinate x = X ∈ [0, 1] to which a crosssectional frame F(X) = (O, t 1 , t 2 , t 3 )(X) is attached, where O(X) is on the midline, t 1 (X) is a unit normal vector perpendicular to the cross-section, and t 2 (X) and t 3 (X) are the unit vectors spanning the cross-sectional plane (see Fig. 3(b)).

The balance of momentum underlying this model can be recast within equation [START_REF] Kim | Soft robotics: a bioinspired evolution in robotics[END_REF], where:

D becomes ∂/∂t, the partial derivative with respect to time t, In the above equations, ρ is the density of the medium, I is the moment of inertia, A is the area of the cross section, r is the position vector of O(X),

q sb becomes [v, ω] T ,
while n and c are the linear and angular cross-sectional stresses along the beam. Fig. 4-(b) shows the application of the method in a soft robot with soft appendices for propulsion in water.

Tendon and thin fluidic actuators can be easily modeled in this framework as active internal wrenches C int that act directly on the cross-sectional stress.

Similar equations can be derived for Cosserat shells [START_REF] Boyer | Poincaré's equations for cosserat media: Application to shells[END_REF]. This modeling framework is not suitable for inflated chambers and some kind of smart material actuations, that inherently require a 3D description of the soft body.

Solution strategies for Cosserat-like models require the use of numerical methods, as they belong to the realm of continuum mechanics. To this end, methods introduced in the previous section are all suitable, i.e., finite and spectral element methods, and they are usually formulated in generalized coordinates. It is worth noticing that Cosserat models, regardless of the resolution strategy, are well suited to hybrid soft-rigid systems since the kinematics of both, rigid bodies and Cosserat material points, take values in SE(3). Furthermore, more recently a serial, strain-based parametrization of the Cosserat rod has been proposed that seamlessly extend the traditional modeling approach of tree-like rigid robots to hybrid soft-rigid system [START_REF] Renda | A geometric and unified approach for modeling soft-rigid multi-body systems with lumped and distributed degrees of freedom[END_REF].

In addition to the many research papers based on the Cosserat approach, a few software and toolboxes have been proposed recently to solve the related equations. Elastica [START_REF] Gazzola | Forward and inverse problems in the mechanics of soft filaments[END_REF] is a free and open-source software, which provides a discrete differential geometry approximation of the Cosserat rod model.

SoRoSim [START_REF] Mathew | Sorosim: a matlab toolbox for soft robotics based on the geometric variable-strain approach[END_REF] is a MATLAB toolbox that implements a strain-parametrization of soft robotic arms. A piecewise-constant strain approximation has also been included as a plug-in in SOFA [START_REF] Adagolodjo | Coupling numerical deformable models in global and reduced coordinates for the simulation of the direct and the inverse kinematics of soft robots[END_REF].

Although Cosserat-like rod and shell models guarantee a good accuracy for soft robots with one or two main dimensions (e.g., slender bodies) undergoing finite deformation, they are not suitable for full 3D soft robot body deformations.

Finite-dimensional parametrization models

The previous models derive from continuum solid mechanics, therefore they lead to sets of partial differential equations, as described in sections 2.2.1 and 2.2.2, that are computationally expensive to solve. It is also possible to describe the behaviour of the soft body via finite dimensional models. These rely on a description of the soft body that leverages a suitably parametrized central axis or "backbone curve", that assumes a prescribed expression. This description leads to sets of ordinary differential equations in contrast to partial differential equations produced by the previous two modeling strategies.

The parametrization of the backbone curve can be defined as follows. Each frame F(X) described in the previous section can be viewed as an element of the group of special Euclidean transformations (i.e., rigid-body displacements).

The origin of this frame traces out the backbone curve as the curvilinear coordinate X varies. The tangent to this curve is the normal to the cross-sectional plane, t 1 (X) (see Fig. 3(c)). For each fixed value of time, the combined rigidbody instantaneous translational and rotational rate of change with respect to arclength is defined by matrix Θ(X) = [F(X)] -1 (∂F(X)/∂X) where F(X) is the roto-translational matrix describing the orientation and translation of the frame (relative to the world frame at the base of the robot), with components R(X) and p(X), respectively [START_REF] Fu | Lateral oscillation and body compliance help snakes and snake robots stably traverse large, smooth obstacles[END_REF][START_REF] Chirikjian | A modal approach to hyper-redundant manipulator kinematics[END_REF][START_REF] Kim | Optimizing curvature sensor placement for fast, accurate shape sensing of continuum robots[END_REF]. In the case when the backbone curve is inextensible (not stretchable), then there is coupling between the rotation matrix and the position vector. For instance, if the tangent to the soft robot arm at its base is e 1 = [1, 0, 0] T , then the tangent at X will be R(X)e 1 and integrating the tangent along the curve generates the position as p(X) = X 0 R(s)e 1 ds . In general, it is always possible to use Θ(X), that is a coordinate-free parametrization. However, when there is coupling between the rotation matrix and the position vector, one can use an alternative parametrization of the backbone curve by expanding rates of rotation parameters such as Euler angles. In the planar case, they reduce to the same parametrization.

The matrix Θ(X) has embedded in it information about instantaneous rotational and translational changes as a function of arclength. These can be extracted to form a vector q sb = [v ω] T , where now v and ω are linear and angular velocity with respect to arclength X instead of time t (as in section 2.2.2) . The vector q sb = [v ω] T can then be expanded using a modal approach, that describes the spatio-temporal behavior of v and ω. This modal expansion is commonly defined as q sb (X, t)

= n i=1 Φ i (X)a i (t), where Φ i (X)
are modes describing the spatial behaviour of the system, and a i (t) are coefficients describing its evolution in time. The v-part of this vector describes how the backbone curve stretches and how adjacent planar sections shear relative to each other. The ω-part of this vector describes bending and twisting, and is related to the classical concepts of curvature and torsion of space curves.

The modal expansion includes as a special case piecewise-constant curvature models by choosing some functions Φ i (X) to be piecewise constant for some values of X and equal to zero for others, as described in [START_REF] Chirikjian | A modal approach to hyper-redundant manipulator kinematics[END_REF].

The parametrization and modal expansion introduced lead to a set of ordinary differential equations, that can be written using the formalism introduced in equation ( 1), where D becomes d/dt, the total derivative with respect to time t, q sb becomes [v ω] T , and

N becomes -1/M [D(s sb , q sb ) + K(s sb )]
, where s sb is the displacement of the soft body.

The nonlinear term just introduced is derived from the robot dynamic model, where M is the inertial matrix, D is a dissipative term that includes internal friction and other dissipative forces (such as Coriolis forces), and K is an elastic term that encapsulate the stiffness of the system.

The coupling with actuation is achieved via an appropriate choice of the modal expansion, and it is commonly represented via C int = α, where α is the vector of actuation forces. When motors specify the angles between rigid links of a traditional robot, the modes would be Dirac delta functions in ω to describe revolute joints or in v to describe prismatic joints, and the weights would be the joint angles, and the actuation forces would be joint torques. In this sense, the modal approach can be used not only for soft/continuum robots but also for ones composed of rigid links, as well as hybrids.

The solution of the models in sections 2.2.1 and 2.2.2 requires the use of numerical methods, e.g., finite and spectral element approximations. In contrast, solution methods for the ordinary differential equations arising in this section are based on iterative schemes that update the weights a i (t) in time. These involve inverse Jacobian iterations as described in [START_REF] Fu | Lateral oscillation and body compliance help snakes and snake robots stably traverse large, smooth obstacles[END_REF][START_REF] Chirikjian | A modal approach to hyper-redundant manipulator kinematics[END_REF][START_REF] Kim | Optimizing curvature sensor placement for fast, accurate shape sensing of continuum robots[END_REF].

The computational benefit of this approach is that instead of solving partial differential equations, it averages properties over each section and reduces the problem to a finite number of ordinary differential equations. This approach has been applied to hyper-redundant manipulators of all kinds (i.e., those consisting of a large number of rigid links, continuum filament, soft) over the past several decades. Successful implementation in the realm of soft robotics applications exist, see for instance [START_REF] Suzumori | Development of flexible microactuator and its applications to robotic mechanisms[END_REF][START_REF] Webster | Design and kinematic modeling of constant curvature continuum robots: A review[END_REF].

Data-driven and machine learning models

An alternative approach to the models above is the use of machine learning techniques. In robotics, learning has been widely used for building the kinematic and dynamic relations between a robot actuators and its position in space. Neural networks have been used to this aim, mainly for control purposes, implementing so-called neuro-controllers. In the same way, in soft robots, neural networks can be used to learn the soft body dynamics and solve the transformations from the actuator space to the space of the soft body deformation and position. Following the formalism in equation ( 1), one tries to learn the nonlinear and the coupling term for internal interactions, N and C int . This results in learning the right-hand side of equation ( 1) in one go, with e.g., a neural network (NN):

Dq sb = NN[q sb (k) → q sb (k + 1)],
where k are different time instances.

In [START_REF] Thuruthel | Control strategies for soft robotic manipulators: A survey[END_REF], an analysis is given of how learning-based blocks can replace some of the steps of the longer transformation chain involved in soft robot control. Diverse network topologies and learning paradigms can be used for this purpose. In [START_REF] Kim | Review of machine learning methods in soft robotics[END_REF], a quantitative survey is presented, showing how supervised learning is more widely used than unsupervised or reinforcement learning, with references to the different techniques used. The survey presented in [START_REF] Wang | A survey for machine learning-based control of continuum robots[END_REF] describes the learning techniques used for obtaining the mapping from actuation space to task space in continuum robots.

An interesting comparison between a model-based and a learning-based approach to the control of a same soft robot arm is presented in [START_REF] Giorelli | Neural network and jacobian method for solving the inverse statics of a cable-driven soft arm with nonconstant curvature[END_REF]. While the model-based method is more accurate in controlling the end-effector position, as resulting in simulation, its error tends to increase together with the model inaccuracies, e.g., fabrication inaccuracies. The learning-based controller error tends to be insensitive to such modeling inaccuracies, especially if training is accomplished on the physical robot arm.

Looking ahead, we envisage a promising direction of progress in the integration of modeling techniques with machine learning, where models can feed neural networks and learning can replace specific sub-system mappings.

More broadly, there is a wealth of literature on data-driven and reducedorder modeling techniques that may be used to capture challenging aspects of soft-robot internal interactions, such as the two-way coupling in fluid structure interactions. System identification [START_REF] Juang | Applied System Identification[END_REF] may be viewed as a predecessor of many modern techniques in machine learning, and this field has a rich history going back at least seven decades [START_REF] Brunton | Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control[END_REF]. For example, considerable work has gone into linear system identification, based on frequency-response and impulse-response data, which is useful for modeling small deflections, especially when feedback control is used to keep the system in a region where linearization is valid.

For large-amplitude motions and unsteadiness, nonlinear model reduction and system identification have become mature techniques [START_REF] Benner | A survey of projection-based model reduction methods for parametric dynamical systems[END_REF][START_REF] Peherstorfer | Data-driven operator inference for nonintrusive projection-based model reduction[END_REF][START_REF] Brunton | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF][START_REF] Loiseau | Constrained sparse Galerkin regression[END_REF][START_REF] Qian | Lift & learn: Physicsinformed machine learning for large-scale nonlinear dynamical systems[END_REF].

538 x y z X X ℱ 0 (X ) = (O 0 , t 1 , t 2 , t 3 )(X ) ℱ s 1 (X ) = (O s 1 , t 1 , t 2 , t 3 )(X ) ℱ s 2 (X ) = (O s 2 , t 1 , t 2 , t 3 )(X )
Backbone Backbone

Rigid cross sections

Computational mesh

ℱ 0 (X ) ℱ p 1 (X ) ℱ p 2 (X ) (a) 3D continuum x y z X X ℱ 0 (X ) = (O 0 , t 1 , t 2 , t 3 )(X ) ℱ s 1 (X ) = (O s 1 , t 1 , t 2 , t 3 )(X ) ℱ s 2 (X ) = (O s 2 , t 1 , t 2 , t 3 )(X )
Backbone Backbone

Rigid cross sections

Computational mesh

ℱ 0 (X ) ℱ p 1 (X ) ℱ p 2 (X ) (b) Shells and rods x y z X X ℱ 0 (X ) = (O 0 , t 1 , t 2 , t 3 )(X ) ℱ s 1 (X ) = (O s 1 , t 1 , t 2 , t 3 )(X ) ℱ s 2 (X ) = (O s 2 , t 1 , t 2 , t 3 )(X )
Backbone Backbone

Rigid cross sections

Computational mesh The 3D continuum solid mechanics model is solved using SOFA [START_REF] Duriez | Soft robot modeling, simulation and control in real-time[END_REF], while the Cosserat model is solved using the methodology described in [START_REF] Armanini | Flagellate underwater robotics at macroscale: Design, modeling, and characterization[END_REF]. In both subfigures (a) and (b), we depict the simulation (top) along with the fabricated (bottom) soft body. In the figure, we also present (c) a simplified pseudo-algorithm for internal interactions in the context of continuum models (i.e., 3D continuum solid mechanics and Cosserat method), noting that the specific details depend on the resolution technique adopted.

ℱ 0 (X ) ℱ p 1 (X ) ℱ p 2 (X ) (c) Parametrization
3 External interactions

Challenges

To benefit from embodied intelligence, we require medium-robot interaction modeling that can grasp the emergence of sensory-motor behaviour from the interplay with the surrounding environment. Modeling these external interactions, and coupling them with the internal interactions discussed in the previous section, closes the gap required to modeling the physics of embodied intelligence.

The task when modeling external interactions is to capture the complex, time-varying forces between the actuated soft-robot structures and the surrounding medium. The dynamics of fluids and deformable solids are generally nonlinear and unsteady, and interactions between the soft robot and the medium may create complex feedback loops and hysteresis effects. The key challenges to modeling external interactions are:

1. the multiphysics nature of the problem, that requires substantial interdisciplinary efforts, and 2. the partially known interfacial physics and unmodeled dynamics (e.g., nonlinear friction and turbulent drag).

In the following, we focus on medium-robot interaction, that we denote with the coupling term C ext introduced in equation ( 1). Here, the surrounding medium is either a fluid or a solid. The robot (i.e. the soft body) is assumed to be described by the models presented in section 2, along with its actuation.

In general, any of the soft-body models in section 2 can be used in conjunction with the models introduced next, unless otherwise specified.

Fluid models

When fluid forces acting on a robot have a tangible impact on movement, as with the octopus in Fig. 2, modeling fluid-robot interaction becomes essential.

Here the task is to capture the complex, time-varying, bidirectional forces between the actuated soft robot structures and the surrounding fluid flow. To this end, the field of fluid-structure interaction [START_REF] Dowell | Modeling of fluid-structure interaction[END_REF] is a key cross-discipline that should be taken into account.

We describe three main modeling strategies, that can be adopted in the context of fluid-robot interaction. These include continuum fluid mechanics models, simplified lumped parameter models, and the relatively new field of machine learning for flow modeling. We remark that we keep the description at a general level, and we avoid entering into details that may hinder grasping the essence of these modeling strategies.

Continuum fluid mechanics models

In analogy to continuum solid mechanics, the equations governing flow physics in a continuum setting can be considered as the most general model of a fluid interacting with the soft body described in section 2 (see Fig. 5(a)).

The equations describing continuum fluid mechanics are the Navier-Stokes equations, constituted, in the most general case, by conservation of mass, momentum and energy. From a soft robot perspective, the action of the fluid modeled by the Navier-Stokes equations emerges in the form of equality constraints at the interface between the soft body and the fluid, denoted by Γ sb,f . In particular, we can write these interface constraints as follows:

I e,f =      q sb = q f on Γ sb,f σ sb • n = σ f • n on Γ sb,f x sb = x f on Γ sb,f . (2a) (2b) (2c)
In equation ( 2), q sb = v sb corresponds to the velocity of the soft body, and q f = v f is the velocity of the fluid at the interface Γ sb,f , n is the normal at any soft body location, σ f is the fluid Cauchy stress tensor composed of the sum of the deviatoric stress tensor τ (f) (that accounts for the viscosity) and a pressure term -pI, and x sb and x f are the positions of the soft body and fluid interfaces, respectively. The solution for the solid velocity is achieved through the modeling strategies introduced in section 2. The solution for the fluid velocity is instead obtained by solving the Navier-Stokes equations. For the purpose of this article, we report the momentum equation, that can be written as follows: the nature of the coupling, the majority of the soft body models presented in section 2 may be used, with the due approximations required for coupling. For instance, in the case of continuum solid mechanics (section 2.2.1), the compatibility conditions are applied to the interface nodes between the soft body mesh and the fluid mesh (or a proxy representing the interface, if the soft body / fluid nodes are not matching). In this case, there exist two mainstream approaches for coupling: partitioned and monolithic. The former uses a structural and a fluid solver iteratively until convergence, while the latter solves a fully coupled system of equations in one go [START_REF] Hou | Numerical methods for fluid-structure interaction-a review[END_REF]. When using continuum fluid mechanics in conjunction with simplified methods, such as finite-dimensional parametrization models (section 2.2.3), one instead needs to lump the effect of the fluid solution into a set of degrees of freedom compatible with the finite expansions used to describe the deformation of the body.

∂(ρ f q f )/∂t + ∇ • (ρ f q f ⊗ q f ) = ∇ • σ f + β f ,
In terms of solution strategies, there exist well-established approaches in the realm of fluid-structure interaction methods. The two main approaches are conforming and non-conforming methods [START_REF] Hou | Numerical methods for fluid-structure interaction-a review[END_REF]. The former enforces equations (2b) and (2c), while the latter enforces equation (2a). Conforming methods can be used in conjunction with partitioned methods, where a popular algorithm is the arbitrary Lagrangian Eulerian approach [START_REF] Souli | Arbitrary Lagrangian Eulerian and Fluidstructure Interaction: Numerical Simulation[END_REF]. Nonconforming methods are typically used with monolithic methods, where a popular algorithm is the immersed boundary method [START_REF] Mittal | Immersed boundary methods[END_REF][START_REF] Taira | The immersed boundary method: a projection approach[END_REF][START_REF] Goza | A strongly-coupled immersed-boundary formulation for thin elastic structures[END_REF].

Continuum fluid mechanics and its interaction with soft bodies have been investigated in biology [START_REF] Dickinson | How animals move: an integrative view[END_REF][START_REF] Lauder | Fish locomotion: recent advances and new directions[END_REF]. Yet, muscular activation and its full coupling with the fluid is still problematic. This issue is also present for soft robots, where muscular activation is replaced by the actuation strategies introduced in section 2.1. For instance, in the case of fluidic actuation, it is necessary to consider both the coupling of the external and of the internal fluid with the soft body. It is exactly this aspect, namely the accurate coupling of the external fluid with the soft body dynamics and actuation system, as well as the coupling of the actuation system with the soft body, that constitute an important modeling challenge in this context. In addition, the demanding computational costs further hamper the use of continuum fluid mechanics models.

To solve the equations arising in this modeling strategy, there exist both commercial and open-source multiphysics packages. These include ANSYS [START_REF]ANSYS software[END_REF] and COMSOL [START_REF]COMSOL software[END_REF], Altair [START_REF]Altair software[END_REF], Flow3D [START_REF]FLOW3D software[END_REF], Adina [START_REF]Adina software[END_REF], and OpenFOAM [START_REF] Weller | A tensorial approach to computational continuum mechanics using object-oriented techniques[END_REF],

to cite a few. However, these are not readily tuned for soft robotics applications.

For example, implementing actuation strategies in a nonlinear elastic material interacting with an external fluid described by the Navier-Stokes equations, where both internal and external interactions are fully coupled with the soft body, remains challenging.

Continuum fluid mechanics can fully capture the effect that a fluid has on a soft robot. However, its inherent computational costs make these models prohibitive to use for practical design and control purposes, and they have found limited use in soft robotics, as of today. As computational resources become more available and underlying algorithms more efficient, the scope of these models will expand and they may be eventually used for design and control

purposes. Yet, today they can already be adopted for high-fidelity simulations that can grasp the complex behaviour of biological and soft robotics systems interacting with a fluid, and extract the key physical principles underlying their embodied intelligence. These principles can in turn improve simplified and less expensive models used for design and control.

Lumped parameter models

In some circumstances, it is acceptable and convenient to simplify the description of the fluid interacting with the soft body. In these cases, one no longer solves the equations governing continuum fluid mechanics, but aggregates the overall effect of the fluid into a set of lumped contributions. The main lumped contributions consist of added mass, drag/lift forces, and buoyancy (see Fig. 5(b)). These can be included as forces into the coupling term C ext as follows: C ext = f added mass + f drag + f lift + f buoyancy , where the subscript of each term is self-explanatory. The added mass, f added mass , is the virtual mass or inertia added to the soft robot due its need to move the fluid surrounding it. For instance, if the soft body was a simple sphere immersed in an incompressible fluid, the added mass would be equal to

f added mass = (1/2)ρ f V sb [Dq f /Dt -dq sb /dt],
where Dq f /Dt is the material derivative of the fluid velocity, dq sb /dt is the total derivative of the spherical soft body velocity, ρ f is the fluid density and V sb is the volume of the spherical soft robot. Drag and lift forces are typically proportional to the velocity of the fluid flowing around the soft body, and they can be calculated as f drag = (1/2)ρ f q f A sb C drag , and

f lift = (1/2)ρ f q f A sb C lift ,
where A sb is the area of the soft body exposed to the fluid, and C drag and C lift are the drag and lift coefficients that, for simplified geometries, are commonly tabulated. Finally, the buoyancy term is typically accounted for as a net vertical force composed by buoyancy plus gravity.

The term C ext just defined enters into the right-hand side of the soft-body equation ( 1) as a forcing term, and can be used in conjunction with any of the models presented in section 2. Obviously, the lumped contributions need to correctly interface with the degrees of freedom of the soft body. In the case of continuum solid mechanics models, (i.e., the ones introduced in sections 2.2.1 and 2.2.2), the lumped forces need to be distributed onto the nodes underlying the discretization of the soft body. Similarly, for finite-dimensional parametrization models (i.e., the ones introduced in section 2.2.3), C ext needs to be applied to the degrees of freedom of the functional parametrization adopted.

The coupling term C ext , and its constitutive components, (i.e., added mass, lift, drag and buoyancy) depend on the shape of the soft body, and therefore they can change over time as effect of actuations or interactions with the environment. Employing such changes to increase efficiency, to direct behaviours, or to improve performances is key to establish quantitative advantages of soft robotics.

In practice, the forces expressed in C ext are often coupled with continuum solid models of soft robots (i.e., the ones described in section 2.2.1 and 2.2.2), such that a complete coupling between continuum fluid mechanics models and continuum solid mechanical models is generally not performed due to its computational costs. Successful implementation of these strategies in soft robotics has been performed, for example in [START_REF] Renda | A unified multi-soft-body dynamic model for underwater soft robots[END_REF] and [START_REF] Armanini | Flagellate underwater robotics at macroscale: Design, modeling, and characterization[END_REF]. In Fig. 6-(a), we

show an example of lumped-parameter model for fluid-solid interaction, used for simulating a flagellate [START_REF] Armanini | Flagellate underwater robotics at macroscale: Design, modeling, and characterization[END_REF].

Eventually, even if external forces are included without modelling the continuum of the fluid, the dependence of such forces on shape-varying coefficients allows reaching a trade-off between accuracy and relevance of simulation, with the computational power available for modelling. Depending on the goal of the simulation, this is similar to fluid-structure interaction in the aerospace industry: while complete coupled simulations are used for design, flight simulators employ a similar coupling to match fidelity of simulation with real-time computation.

Data-driven and machine learning models

Recent data-driven and machine learning methods have become widely used for modeling complex fluids [START_REF] Duraisamy | Turbulence modeling in the age of data[END_REF][START_REF] Brunton | Machine learning for fluid mechanics[END_REF][START_REF] Vinuesa | The potential of machine learning to enhance computational fluid dynamics[END_REF]. Several recent approaches have leveraged machine learning to direct speed up high-fidelity simulation of the Navier-Stokes equations, especially the ones involving turbulence [START_REF] Bar-Sinai | Learning datadriven discretizations for partial differential equations[END_REF][START_REF] Kochkov | Machine learning accelerated computational fluid dynamics[END_REF][START_REF] Wang | Incorporating symmetry into deep dynamics models for improved generalization[END_REF][START_REF] Li | Fourier neural operator for parametric partial differential equations[END_REF]. Indeed, for complex flows, it is often impractical to resolve all scales of the flow, and instead researchers employ turbulence models, such as the Reynolds averaged Navier-Stokes (RANS) equations or large eddy simulation (LES). These approximate the smaller scales, and allow for less computationally expensive simulations of the Navier-Stokes equations. Machine learning is rapidly advancing these computational fields [START_REF] Duraisamy | Turbulence modeling in the age of data[END_REF], providing enhanced data-driven turbulence closures [START_REF] Ling | Reynolds averaged turbulence modelling using deep neural networks with embedded invariance[END_REF][START_REF] Maulik | Subgrid modelling for two-dimensional turbulence using neural networks[END_REF][START_REF] Novati | Automating turbulence modelling by multi-agent reinforcement learning[END_REF][START_REF] Beetham | Formulating turbulence closures using sparse regression with embedded form invariance[END_REF][START_REF] Beetham | Sparse identification of multiphase turbulence closures for coupled fluid-particle flows[END_REF]. For even further reduction in computation, it is often possible to develop reduced-order models (ROMs) that are tailored to a specific flow and provide an optimal balance between accuracy and efficiency. Reducedorder models are typically at least partially data-driven, as they are based on modal decompositions [START_REF] Taira | Modal analysis of fluid flows: An overview[END_REF], and they have close connections to machine learning. Several recent approaches have provided more accurate and generalizable ROMs, for example based on sparse regression [START_REF] Brunton | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF][START_REF] Loiseau | Constrained sparse Galerkin regression[END_REF][START_REF] Loiseau | Sparse reduced-order modeling: sensor-based dynamics to full-state estimation[END_REF][START_REF] Deng | Low-order model for successive bifurcations of the fluidic pinball[END_REF][START_REF] Deng | Galerkin force model for transient and post-transient dynamics of the fluidic pinball[END_REF], operator inference [START_REF] Peherstorfer | Data-driven operator inference for nonintrusive projection-based model reduction[END_REF][START_REF] Qian | Lift & learn: Physicsinformed machine learning for large-scale nonlinear dynamical systems[END_REF], and the use of deep neural networks to learn effective coordinate systems [START_REF] Lee | Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders[END_REF]. Other promising recent techniques for physics computation based on machine learning include physics-informed neural networks [START_REF] Raissi | Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations[END_REF] and deep operator networks [START_REF] Lu | Learning nonlinear operators via deeponet based on the universal approximation theorem of operators[END_REF].

In addition to using machine learning techniques to develop surrogate models for the complex fluid environment, there is also a need to model discrepancies that arise in complex multiphysics applications. The equations of a simple fluid are relatively well understood, but interfacial dynamics, non-Newtonian fluids, and multiphase flows all pose significant modeling challenges even to represent the fundamental physical effects in a full-fidelity model, let alone in a reduced-order model. There are considerable efforts to use machine learning approaches to model these discrepancies between observed data and idealized physical models. More generally, machine learning is complementing the mature fields of system identification and reduced-order modelling [START_REF] Benner | A survey of projection-based model reduction methods for parametric dynamical systems[END_REF][START_REF] Peherstorfer | Data-driven operator inference for nonintrusive projection-based model reduction[END_REF][START_REF] Brunton | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF][START_REF] Loiseau | Constrained sparse Galerkin regression[END_REF][START_REF] Qian | Lift & learn: Physicsinformed machine learning for large-scale nonlinear dynamical systems[END_REF],

where there are a host of modeling techniques ranging from rather opaque (i.e., black box) to transparent (i.e., white box). There have been several recent efforts to augment physics-based approaches with machine learning to model the discrepancy between the physics model and observed data; this is often referred to as residual learning and has been used to model fluid effects in flying robots [START_REF] Shi | Neural lander: Stable drone landing control using learned dynamics[END_REF].

Solid models

A model of solid-robot interaction should accurately capture how normal forces, tangential forces or friction, and adhesion/cohesion forces affect and can therefore be leveraged by the soft body. To this end, contact mechanics [START_REF] Johnson | Contact Mechanics[END_REF] and tribology [START_REF] Vakis | Modeling and simulation in tribology across scales: An overview[END_REF] represent the key broad cross-disciplines for modeling solid-robot interaction, where modeling the microstructure, e.g., roughness and topography, can be of crucial importance for tasks like gripping [START_REF] Dalvi | Linking energy loss in soft adhesion to surface roughness[END_REF].

In the following, we outline the main modeling strategies that can be adopted in the context of solid-robot interaction. These include the use of continuum solid mechanics to model a deformable solid, simplified lumped parameter models, and more recent machine learning strategies. A key for these models to be accurate is related to the frictional interaction, and consequent constraints at the interface between the solid and the soft robot.

Continuum three-dimensional solid mechanics models

The models underlying continuum solid mechanics have been introduced in section 2.2.1, albeit for the soft robot. These models can be used also to describe a surrounding deformable solid medium (see Fig. 5(c)). Its interaction with the soft robot is described by the coupling term C ext introduced in equation ( 1). This is typically the result of a series of constraints that need to be satisfied at the interface between the soft robot and the solid medium. Such constraints are similar to the one introduced in section 3.2.1, except that they need to account for complementarities, that is, when the system transition from no contact to contact, or from static friction (i.e., no tangential relative movement between two elastic surfaces) to dynamic fiction (i.e., relative tangential movement).

For instance, the normal contact between two elastic surfaces can be modeled using the Signorini's framework, that uses a distance function S n to measure the distance between two surfaces along the normal direction. Signorini's formulation can be written as S n (q) ⊥ σ n , where σ n is the stress along the normal direction n at the contact interface, and the symbol ⊥ denotes complementarity: if S n (q) > 0 then σ n = 0 (no contact); if S n (q) = 0 then σ n > 0 (contact). Tangential or frictional contact can instead be modeled using the Coulomb's law framework, where the complementarity is given by the relation

σ t = µ s σ n .
Here, complementarity arises from static vs. dynamic friction. For σ t ≤ σ max there is static friction (no relative sliding of the two surfaces), while for σ t > σ max there is dynamic friction (relative sliding of the two surfaces), where σ max is the maximum value of stress that allow static friction.

These interface conditions need to be encapsulated into the continuum solid mechanics models describing the soft robot and the elastic solid medium interacting with the robot, in a similar manner as described in section 3.2.1, for the constraints in equation [START_REF] Pfeifer | Self-organization, embodiment, and biologically inspired robotics[END_REF]. They can for instance be embedded into the equations through the term C ext = I e,s Λ, using Lagrange multipliers for the constraints. However, in contrast to section 3.2.1, one needs to account for the complementarities in the contact and friction laws. These complementarities create jumps in velocities. When the coupled system transitions from one regime to another (e.g., from no contact to contact, or from static to dynamic friction), the inertial terms underlying the system of equations is not defined.

In order to solve this contact problem, it is necessary to rely on numerical methods that account for these complementarities, and the singularities they produce. To this end, a strategy is first to calculate S n (q) using a measurement of proximity distance, interpenetration distance, interpenetration volumes, or a precise measurement of the moment and the contact configuration between two simulation steps. Following this detection, a set of constraints are defined, often on contact points (but one can extend to volumes or to non-planar surfaces).

Finally, the solution can be obtained with different numerical strategies such as Lagrange multipliers, penalty methods and augmented Lagrangian methods, in conjunction with optimization solvers dedicated to complementarity problems.

The time-integration of the resulting equations is achieved via event-driven or time-stepping methods that account for the singularity in the underlying equations -see e.g., [START_REF] Studer | Simulation of non-smooth mechanical systems with many unilateral constraints[END_REF], [START_REF] Acary | Projected event-capturing time-stepping schemes for nonsmooth mechanical systems with unilateral contact and coulomb's friction[END_REF].

An example of successful implementation of this strategy can be found in [START_REF] Coevoet | Optimization-based inverse model of soft robots with contact handling[END_REF]. Here, the contact problem is formulated to solve an inverse model for the control of a soft robot, which leads to writing a quadratic problem with complementarity constraints (also referred to as QPCC).

Most of the major multiphysics computational software previously cited, including ABAQUS [START_REF]ABAQUS software[END_REF], ANSYS [START_REF]ANSYS software[END_REF] and COMSOL [START_REF]COMSOL software[END_REF], can be used to solve contact between two elastic solids, yet they are not necessarily tailored to soft robotics applications. In the article [START_REF] Collins | A review of physics simulators for robotic applications[END_REF], a literature review of physical engines for simulation in robotic applications is listed, most of these engines support contact modeling. On the other hand, not all engines are adapted to deformable robots, in particular to simulate volume deformations. We note the SOFA software [START_REF]sofa framework[END_REF], originally dedicated to medical simulation, offers plug-ins dedicated to soft robotics and proposes implementations of FEM and Cosserat rods that are compatible with contact modeling. An example of these capabilities is depicted in Fig. 6-(b), where SOFA simulates a finger touching a deformable solid.

Similar considerations as for continuum fluid mechanics can be made. In particular, these models are computationally expensive, and they may lead to impractical solutions for soft robotics applications today, especially in the context of real-time control.

Lumped parameter models

Similarly to what we have introduced for fluid models, we can consider a simplified approach, where we can aggregate the effect of contact between two surfaces into lumped contributions. For instance, when the external solid can be considered rigid, one can introduce some simplifications to the model outlined in section 3.3.1. In particular, the state q s of a rigid body can be represented by the position of its center of mass, its orientation, and its linear and angular velocity (see Fig. 5(d)). Once these are known, it is possible to formulate the rigid-robot interaction problem. This implies identifying the point of contact between the soft body and the rigid solid, and applying an equivalent frictional force f and torque τ , induced by the interaction between the soft and the rigid body. The interaction surface between the rigid solid and the soft robot is typically modeled as a planar surface. This leads to models that contain only the three degrees of freedom associated to frictional forces at the planar surface, f .

In practice, a soft robot interacting with a rigid body leads to non-planar contact surfaces, with multiple points of contact. This introduces three additional degrees of freedom, and provide a six-dimensional model for the normal and frictional wrenches, as both force f and torque τ are three-dimensional.

Following the formalism adopted in [START_REF] Xu | 6dls: Modeling nonplanar frictional surface contacts for grasping using 6-d limit surfaces[END_REF], one can define the normal force and torque as f n = -S p • ndS, and f n = -S p • [(r × n]dS, respectively, and the frictional force and torque as f t = -µ S p • v r dS, and

τ t = -µ S p•[r × v r ]
dS, respectively, where S is the contact surface, v r is the relative velocity between the rigid solid and the soft robot, n is the normal to the contact surface, r is the torque arm, p is the contact pressure distribution, and µ is the friction coefficient.

The coupling is achieved by encapsulating these forces and torques into the coupling term C ext . This acts as a forcing term on the right-hand side of the soft body equation, that is

C ext = f n + τ n + f t + τ t .
The computational costs associated with these models are relatively small compared to the interaction between two deformable solids.

Applications of this modeling strategies in soft robotics exists, both in terms of analytical approaches [START_REF] Xydas | Modeling of contact mechanics and friction limit surfaces for soft fingers in robotics, with experimental results[END_REF][START_REF] Majidi | Influence of surface traction on soft robot undulation[END_REF] and computational models [START_REF] Duriez | Soft robot modeling, simulation and control in real-time[END_REF][START_REF] Todorov | Mujoco: A physics engine for model-based control[END_REF][START_REF] Deimel | A novel type of compliant and underactuated robotic hand for dexterous grasping[END_REF],

that focused on simulating contact behaviours of soft robot grasping and manipulation as well as crawling [START_REF] Lipson | Automatic design and manufacture of robotic lifeforms[END_REF][START_REF] Hiller | Automatic design and manufacture of soft robots[END_REF]. Inverse model of soft robot in situation of frictionless contact [START_REF] Coevoet | Optimization-based inverse model of soft robots with contact handling[END_REF] and adhesive contact situations has been investigated for control of manipulation and locomotion [START_REF] Coevoet | Soft robots locomotion and manipulation control using fem simulation and quadratic programming[END_REF].

These models can be broadly applied for describing several tasks involving the interaction of the soft robot with e.g., rigid solids, and similar considerations as in section 3.2.2 can be made.

Data-driven and machine learning models

Frictional and normal contact modeling is notoriously difficult, and a high prediction accuracy is onerous to achieve, especially for soft-bodied robots that interact with environments or objects with vastly different surface properties.

Within this context machine learning and data-driven approaches have found limited use thus far, with few exceptions.

Recent efforts explore the augmentation of simulation representations with neural networks to close sim-to-real gaps (see, e.g., [START_REF] Hwangbo | Learning agile and dynamic motor skills for legged robots[END_REF][START_REF] Golemo | Sim-to-real transfer with neural-augmented robot simulation[END_REF][START_REF] Battaglia | Interaction networks for learning about objects, relations and physics[END_REF][START_REF] Jiang | Data-augmented contact model for rigid body simulation[END_REF][START_REF] Sundaram | Learning the signatures of the human grasp using a scalable tactile glove[END_REF]). These efforts attempt to reduce the uncertainty that contact with unknown objects or environments introduces. A naive approach is to learn the mapping of the previous to the next state directly. This would mean that the nonlinear equations N in equation [START_REF] Kim | Soft robotics: a bioinspired evolution in robotics[END_REF], and also the actuation and solid-solid interaction constraint force terms, are replaced with a neural network.

However, this is all but practical for complex robots, and it is best to acquire data for interactions, such as frictional contact, that lead to prediction inaccuracies in our state-of-the-art models, and to only learn corrections.

More specifically, we could learn corrections for the frictional forces that act on the soft-body, improving a first-order estimate that we get from a known simulation representation.

Continuum fluid mechanics (Navier-Stokes equations) ) solid-solid interaction using 3D continuum mechanics. The details of the fluid-solid interaction simulation are described in [START_REF] Armanini | Flagellate underwater robotics at macroscale: Design, modeling, and characterization[END_REF], while the solid-solid interaction simulation has been solved using SOFA [START_REF] Duriez | Soft robot modeling, simulation and control in real-time[END_REF]. In both subfigures (a) and (b), we depict the simulation along with the fabricated soft body. In the figure, we also present (c) a simplified pseudo-algorithm for fluid-solid and solid-solid interactions in the context of continuum mechanics models (i.e., 3D continuum fluid and solid mechanics), noting that the specific details depend on the resolution technique adopted.

∂ ( ρ f q f) /∂t + ∇ ⋅ ( ρ f q f ⊗ q f) = ∇ ⋅ σ f + β f Compatibility constraints at interface Γ sb,f ℐ e,f = q sb = q f σ sb ⋅ n = σ f ⋅ n x sb = x f ext = Λ T ℐ e,f Γ sb,f (
sb = q s σ sb ⋅ n = σ s ⋅ n x sb = x s ext = Λ T ℐ e,s S n (q) ⊥ σ n σ t = μ σ n Γ sb,s ∂q s /∂t = ∇ ⋅ σ s σ s = σ T s (c)

From models to practice

As emerging from the previous sections, presenting the models involved in describing soft robots and embodied intelligence inside a unified framework is challenging, yet helpful. Therefore, it is not surprising that a structured modeling framework for the physics of interactions in soft robotics remains elusive.

The models presented range from high-dimensional (sections 2. control and design purposes. The latter have been instead used with different degrees of success due to their simplifications, that may be unable to capture the rich physics of internal and external interactions.

High-dimensional models can provide a framework for simulation of internal and external interactions. Embodied intelligence shifts the challenges of realizing a specific task from the control algorithm to such interactions, and therefore to the design stage. The models presented could be employed to investigate the optimal design to achieve a specific task with minimal feedback control strategies, or open-loop ones. As an example, in case of manipulation, diverse designs could be investigated to obtain effective grasping on a variety of objects by employing the same basic low-level control strategy (by inflating a pneumatic chamber, or pulling a tendon). It is worth to mention that embodied intelligence points toward the effective expression of a certain behavior, and not to the accurate expression of it. That is, the fingers of a soft gripper are free to adapt to external interactions up to the point that the grasp is still stable, and no control actions are expected to place a specific finger in a specific position. Although most of the time this characteristic is obtained by physical prototyping and testing a mock-up ("trial-and-error" approach), it can also be obtained with the high-dimensional models presented, with the further benefit of including stochastic disturbances to the actuation, materials or interaction. The embodied intelligence principles allow the design of robust agents that express the desired behavior in presence of disturbances, and therefore the sim-to-real gap can be reduced by looking for such robust designs.

These simulations, in turn, can be used to understand the principles behind an observed and/or desired behaviour and may eventually be adopted for automatic design [START_REF] Lipson | Automatic design and manufacture of robotic lifeforms[END_REF][START_REF] Hiller | Automatic design and manufacture of soft robots[END_REF][START_REF] Lipson | Challenges and opportunities for design, simulation, and fabrication of soft robots[END_REF] in the future, if computational power and algorithms will allow for a sufficiently fast model-driven workflow. Yet, having these models permits constructing improved approximated and low-dimensional representations of soft robots, that capture their overall behaviour at a cheaper computational cost, feasible for design and control purposes today. Indeed, the high-dimensional models can capture the key components underpinning embodied intelligence physical principles, thereby allowing for a more accurate low-dimensional description of how soft robots can achieve the embodied intelligence of biological systems. As an example, Fig. 2 depicts a low-dimensional model for octopus' underwater locomotion [START_REF] Calisti | Morphological and control criteria for self-stable underwater hopping[END_REF]. A low-dimensional description is typically constituted by a representation of the behaviour of the system in a reduced space, composed by much fewer degrees of freedom than the original high-dimensional counterpart [START_REF] Goury | Fast, generic, and reliable control and simulation of soft robots using model order reduction[END_REF][START_REF] Chenevier | Reduced-order modeling of soft robots[END_REF]. Therefore, one must identify a suitable low-dimensional representation that retains a sufficient accuracy to describe the whole system.

The discovery of low-dimensional representation (often referred to as fundamental or parsimonious models) is currently done empirically by observing the system. Whether a high-dimensional model is available, the synthesis to the lower-level counterpart (if any) can be obtained more systematically. Accurately informed fundamental models can capture the overall behaviour of the soft robot and expose a few key aggregate physical parameters to describe a specific task (e.g., locomotion, arm movements, etc.). These fundamental models provide a generic base upon which specific morphology/control could be developed, following the template/anchor approach proposed in [START_REF] Full | Templates and anchors: neuromechanical hypotheses of legged locomotion on land[END_REF]. The analysis of the fundamental models might increase the comprehension of the system and expose quantitative advantages of compliant robots over rigid ones.

This concept was successfully employed in underwater legged robots, where shape-dependent forces, elastic leg elements, and pushing-based actuation interweave to shape the basin of attraction of the hopping limit cycle. Shape morphing proved how deformable bodies can be exploited to overcome actuation limits [START_REF] Picardi | Bioinspired underwater legged robot for seabed exploration with low environmental disturbance[END_REF]. Another remarkable example is presented in [105], where the elastic components of a soft robotic squid have been exploited to increase swimming efficiency. A second-order forced oscillation model was developed to catch the relationship between internal and external forces, mediated by actuation frequency, in an approximate and simplified way. In both cases, they are not a low-level counterpart of the models presented in section 2 and 3, but they are fundamental models which capture the specific behavior of interest, and can provide insights into the design of actual soft robots.

All three of the main types of internal and external interactions with fluid and solid media may benefit from recent advances in machine learning for modeling dynamical systems [START_REF] Brunton | Machine learning for fluid mechanics[END_REF][START_REF] Bar-Sinai | Learning datadriven discretizations for partial differential equations[END_REF][START_REF] Kochkov | Machine learning accelerated computational fluid dynamics[END_REF]. However, it is important to understand the considerable open challenges associated with using machine learning models for design, optimization, and control. Many machine learning algorithms, such as deep learning, suffer from limited interpretability and generalizability.

Although they can capture patterns in training data, the resulting models are often inscrutable and may have limited ability to model new scenarios that are not captured in the training data. For applications in image classification and natural language translation, training data is extensive and many plausible uses of a machine learning model may be considered as interpolatory in nature. In contrast, when using machine learning to accelerate the design of a new physical system, such as a soft robot, the model must be able to generalize, or extrapolate, to new regimes that have not been characterized in the training data, otherwise a model would not be needed. This capability is particularly important for inverse-based design optimization, where it is necessary to model the dependence of the system on the parameters it will be optimized over. Further, when using machine learning models for design optimization, any computational savings must be balanced against the cost of acquiring training data, which is often significant. Using machine learning for control is easier to justify, as model training is considered an off-line expense that is justified by the increased speed and performance of the model in on-line use. It is important to note that these challenges have been debated for decades, long before the rise of machine learning, to understand the role of reduced-order models more generally in design optimization and control.

Conclusions

We walked through the quest for modelling soft robots with a focus on how to describe the physics involved in their embodied intelligence. We argue that modelling the deformations of a soft robot body under internal actuation forces and external interaction forces can capture the practical essence of embodied intelligence. We outlined the mathematical models used for describing such

Bespoke tools written in various languages

Table 2: Summary of models for external interactions, with reference to equation [START_REF] Kim | Soft robotics: a bioinspired evolution in robotics[END_REF]. As discussed in sections 3.2.3 and 3.3.3, machine learning techniques can be also proposed to solve the same modeling problems.

  the velocity twist field of the soft-body composed of the linear and angular velocity, respectively, N becomes [N sb,v , N sb,ω ], with N sb,v = 1/(ρA)(∂n/∂x + n), and N sb,ω = (I -1 /ρ)(-ω × (ρIω) + ∂c/∂X + ∂r/∂X × n + c).

Fig. 3 :

 3 Fig.3: A conceptual illustration of the models for internal interactions described in section 2, applied to the octopus arm of Fig.2. (a) A representation of the mesh description in 3D continuum mechanics models. (b) Cosserat's approach for rods. (c) Finite-dimensional parametrization models.

Fig. 4 :

 4 Fig.4: Examples of practical simulations in the context of internal interactions achieved via (a) 3D continuum solid mechanics, and (b) Cosserat. The 3D continuum solid mechanics model is solved using SOFA[START_REF] Duriez | Soft robot modeling, simulation and control in real-time[END_REF], while the Cosserat model is solved using the methodology described in[START_REF] Armanini | Flagellate underwater robotics at macroscale: Design, modeling, and characterization[END_REF]. In both subfigures (a) and (b), we depict the simulation (top) along with the fabricated (bottom) soft body. In the figure, we also present (c) a simplified pseudo-algorithm for internal interactions in the context of continuum models (i.e., 3D continuum solid mechanics and Cosserat method), noting that the specific details depend on the resolution technique adopted.

  where ρ f is the fluid density, and β f is a forcing term acting on the fluid (e.g., gravity). Depending on the nature of the problem, we can have several approximations, including incompressible Navier-Stokes equations, where the density is a constant, compressible Navier-Stokes equations, where the density is not constant and the flow might develop shock waves, thin-layer Navier-Stokes equations, that describes flow in thin layers, and potential flow equations that describe irrotational flows, to cite a few. The compatibility conditions in equation (2) can enter into the equations describing the soft body (1) as Lagrange multipliers C ext = Λ T I e,f , or being imposed as external equality constraints I e,f or as boundary conditions. Given

Fig. 5 :

 5 Fig. 5: A conceptual illustration of the models for external interactions described in section 3, applied to the octopus arm of Fig. 2. (a) Continuum fluid mechanics models. (b) Lumped fluid models. (c) Continuum solid mechanics models. (d) Lumped solid models.

Fig. 6 :

 6 Fig.6: Examples of practical simulations in the context of external interactions for (a) fluid-solid interaction using a lumped-parameter model and for (b) solid-solid interaction using 3D continuum mechanics. The details of the fluid-solid interaction simulation are described in[START_REF] Armanini | Flagellate underwater robotics at macroscale: Design, modeling, and characterization[END_REF], while the solid-solid interaction simulation has been solved using SOFA[START_REF] Duriez | Soft robot modeling, simulation and control in real-time[END_REF]. In both subfigures (a) and (b), we depict the simulation along with the fabricated soft body. In the figure, we also present (c) a simplified pseudo-algorithm for fluid-solid and solid-solid interactions in the context of continuum mechanics models (i.e., 3D continuum fluid and solid mechanics), noting that the specific details depend on the resolution technique adopted.

q 2 Fluidq

 2 sb = q f σ sb • n = σ f • n x sb = x f Cext = Λ T I e,f(e.g., Lagrange multipliers) forces applied on the right-hand side of equation (1) as aggregated contributionsf added mass + f drag + f lift + f buoyancy sb = qs σ sb • n = σs • n x sb = xs Cext = Λ T Ie,s (e.g., Lagrange multipliers) Complementarities Sn(q) ⊥ σn σt = µσn SOFA [8] ABAQUS [13] ANSYS [14] COMSOL [15] Lumped parameters Γsb,s section 3.3.2 Normal and tangential forces and torques applied on the right-hand side of equation (1) as aggregated contributions f n , +τ n + f t + τ t SOFA [8] and other bespoke tools written in various languages

  Continuum solid

	Γ sb,s ext = f f n Compatibility constraints at the interface Γ sb,s τ n f t τ t
	Lumped parameter models
	(Solid mechanics)
	(d) Lumped solid

n + τ n + f t + τ t

Table 1 :

 1 Summary of models for soft body mechanics and internal interactions, with reference to equation[START_REF] Kim | Soft robotics: a bioinspired evolution in robotics[END_REF]. As discussed in Section 2.2.4, machine learning techniques can be also proposed to solve the same modeling problems.

		Model	D	q sb	N sb	C int	Software
							tools
		Continuum
		3D solid
		mechanics
		Computational mesh		X		X
		z
		x	y	ℱs 2 (X ) = (Os 2 , t1, t2, t3)(X ) ∂/∂t v sb Backbone	Backbone ∇ • σ sb ℱp 2 (X )	Actuation	SOFA [8]
					ℱs 1 (X ) = (Os 1 , t1, t2, t3)(X )	ℱp 1 (X )	imposed as	ChainQueen [10]
				Rigid cross sections	ℱ0(X ) = (O0, t1, t2, t3)(X )	ℱ0(X )	equality and inequality	Evosoro [12] Voxelyze [11]
		section 2.2.1				constraints: Ie, I i , through Λ T Ie + Θ T I i	ABAQUS [13] ANSYS [14] COMSOL [15] Altair [16]
							(e.g., Lagrange
							multipliers)
							[used with
							1D, 2D, 3D
							actuations]
							[N sb,v , N sb,ω ] T	Actuation imposed as
							N sb,v =	active
							1/(ρA)(∂n/∂x + n)	internal	SOFA [8]
		Rods and shells, e.g., Cosserat X	∂/∂t [v, ω] T	N sb,ω = +∂r/∂X × n + c] +∂c/∂X+ I -1 /ρ[-ω × (ρIω)+	wrenches actuation] 1D, 2D [used with α, τ	Elastica [28] SoRoSim [29]
		Backbone
			ℱs 2 (X )
			ℱs 1 (X )
		Rigid
		cross
		sections
		ℱ0(X )
		section 2.2.2
							-1/M [D(s sb , q sb )+
							+K(s sb )]
		Finite-dimensional parametriza-tion	d/dt [v ω] T	s sb is the displacement of the soft body	Actuator forces grouped in a vector α
	X	X					that acts as
							a forcing
	one	Backbone	ℱp 2 (X )				term on the
							right-hand side
			ℱp 1 (X )				[used with 1D
							actuation]
		ℱ0(X )
		section 2.2.3

(O0, t1, t2, t3)(X ) ℱs 1 (X ) = (Os 1 , t1, t2, t3)(X ) ℱs 2 (X ) = (Os 2 , t1, t2, t3)(X )

We note that the surrounding environment may be constituted of heterogeneous granular material. In this case, one may need to rely on different methods than the ones presented in this paper. Yet, the case of solid or fluid media is sufficiently general to provide an effective guide on modeling interactions between soft bodies and the environment.

Acknowledgments Gianmarco Mengaldo acknowledges NUS support through his start-up grant (R-265-000-A36-133). Gregory Chirikjian acknowledges MOE Tier 1 grant (R-265-000-655-114) . Cecilia Laschi acknowledges NUS support through her start-up grant (R-265-000-A31-133 and R-265-000-A31-731) .Part of this research is supported by the National Research Foundation, Singapore, under its Medium Sized Centre Programme -Centre for Advanced Robotics Technology Innovation (CARTIN). Gianmarco Mengaldo, Gregory Chirikjian, and Cecilia Laschi acknowledge MOE Tier 2 grant 'REBOT'. This work was supported in part by the Khalifa University of Science and Technology under Grants CIRA-2020-074, RC1-2018-KUCARS, and the Office of Naval Research Global under Grant N62909-21-1-2033.

internal and external interactions in a soft robot body, including related software tools, and we discussed how to use them in practice. For the first time, we show a unified view of the multiphysics of interactions arising from embodied intelligence, and we link them to the design of soft robots.

Tables 1, and 2 summarize this analysis (not including data-driven techniques), and provide a practical and concise guide to the modelling methods that can be beneficially used in soft robotics , enabling soft robots to fully leverage on embodied intelligence, acquire unprecedented abilities and respond to unmet needs, ultimately contributing to further soft robotics progress.

We argue that, in contrast to the current trial (physically build the robot) and error (robot testing) approach, soft-robot design can transition to a model-informed workflow. Indeed, prototype-driven design (also referred to as trial-and-error) was likely the fastest and most efficient way to proceed, especially in the soft robotics exploratory phase. Today, we are at a stage where computational modelling within a model-driven umbrella can enable (i) scaling up soft-robot design in response to application needs, (ii) holistic model-based control embedding external interactions, and (iii) high-fidelity simulations, opening the path towards soft-robot digital twins. This transition is within grasp, in an interdisciplinary dialogue of roboticists with communities like computational physics, applied mathematics, scientific computing and machine learning. This unified approach will allow soft robotics to thrive in the next few decades and establish itself as a model-driven scientific discipline that can tangibly impact human activities.