Inverse dynamics of underactuated flexible mechanical systems

An alternative approach to the inverse dynamics of flexible mechanical systems is presented. In contrast to a sequential discretization in space and time a simultaneous space‐time discretization is applied to the problem of inverse dynamics of underactuated systems. In particular, a space‐time finite element formulation will be compared with a formulation which is based on the method of characteristics. Numerical examples are presented which underline the importance of solving this class of problems in space and time simultaneously.


Introduction
Servo constraints have been successfully applied to solve the inverse dynamics of discrete mechanical systems. In this approach the equations governing the motion of the discrete mechanical system at hand are supplemented by algebraic servo constraints. The servo constraints serve the purpose of partially prescribing the motion of the mechanical system (cf. [1,5]). In principle the same approach can also be applied to the inverse dynamics of flexible mechanical systems such as elastic ropes and beams. To this end a discretization in space needs be applied first to generate the discrete mechanical system. Then servo constraints can be appended leading again to differential algebraic equations (DAEs). However, the index of the resulting DAEs can be quite large hindering their numerical solution. In this contribution, an alternative approach will be presented. The present work focuses on mechanical systems whose motion is governed by quasilinear hyperbolic partial differential equations (PDEs). In particular a system is considered whose motion is governed by: Here, s ∈ S = [0, l] contains the arc-length of a reference curve in R n dim , n dim ∈ {1, 2, 3}, t ∈ T = [0, t e ] is the time domain of interest, and r(s, t) ∈ R n dim . We further introduce the space-time domain Ω = S × T . The main task is to find f (t) ∈ R n dim for g(t) ∈ R n dim prescribed, such that (1) is satisfied. In section 2 and 3 a formulation based on the method of characteristics and a space-time finite element formulation is introduced respectively. Here, n dim = 1 and B(s, t) = b 2 are

Method of characteristics
The PDE in (1) can be recasted in the form of a system of first order PDEs by introducing q(s, t) = ∂ t r(s, t) and p(s, t) = ∂ s r(s, t): Using the column vectors x = q p and C = c 0 , alongside the square matrix A ∈ R 2×2 , the system of first order PDEs (2) then reads: In the method of characteristics (cf. [3,4]) a curve t = k(s) is called a characteristic curve if Since the considered problem is hyperbolic, (4) leads to two real-valued solutions for the characteristic curves along which q(s, t) and p(s, t) satisfy the following system of ODEs: The expression ( · ) α , α = 1, 2, refers to the two characteristic curves. The resulting system of ODEs (5), can be solved numerically, e.g. by using finite differences. The boundary and initial conditions specified in (1) can be applied directly at the nodes of the characteristic net.  (1) can be rewritten as a system of first order in time PDEs: Multiplying each equation in (6) with a sufficiently smooth testfunction w 1 (s, t) and w 2 (s, t), respectively, and integrating over the space-time domain Ω, leads together with the Neumann boundary conditions ∂ s r(0, t) = f (t) and ∂ s r(l, t) = 0 after integrating by parts to The servo-constraint r(l, t) = g(t) which has to be satisfied for all t ∈ T can be enforced in a weak sense by using a test function w 3 (t) to obtain T w 3 (t) · (r(l, t) − g(t)) dt = 0.
Equations (3) and (8) constitute the newly proposed space-time finite element formualtion. The test functions w 1 (s, t), w 2 (s, t), w 3 (t), along with the trial functions r(s, t), v(s, t), f (t) can now be approximated by piecewise continuous polynomials (e.g. Lagrangian shape functions). Applying standard finite element procedures yields an algebraic system of equations for the determination of the nodal degrees of freedom associated with the discrete trial functions.

Numerical examples
Linear elastic bar: A bar with length l, cross-sectional area A, density ρ and Young's modulus E is investigated. The task is to find the force F (t) which is acting on one end of the bar (s = 0) such that the other end (s = l) tracks a prescribed trajectory g(t). Assuming linear constitutive relations and linear kinematics, the servo-constrained longitudinal wave propagation in the bar is governed by problem (1), by setting n dim = 1, B(s, t) = E/ρ, c = 0 and f (t) = F (t)/EA. The methods presented in sections 2 and 3 yield numerical results which coincide very well with the analytical reference solution (cf. Fig.1). Here g(t) claims a sinusoidal rest-to-rest maneuver of the bar at s = l. Both methods yield numerical results which coincide very well with the analytical reference solution. For comparison, a spatial discretization of the linear elastic bar using finite elements leads to a semi-discrete servo-constraint problem. The semi-discrete problem can be viewed as spring-mass system in which the number of masses, say n, corresponds to the number of nodes in the finite element model. It can be shown that the index of the DAEs is 2n + 1, see [2,6]. Accordingly, using a reasonably accurate finite element discretization (in space) yields DAEs with excessively high index.
Nonlinear elastic string: The second example deals with large planar (n dim = 2) deformations of an elastic string. In the undeformed stress-free reference configuration the string has length l, cross-sectional area A, density ρ and Young's modulus E. The task is to find the external force F (t) which is acting on one end of the string (s = 0) such that the other end (s = l) follows a prescribed trajectory g(t) ∈ R 2 . The motion of the string is characterized by r(s, t) ∈ R 2 , and the force in the extensible string is denoted by n(s, t) ∈ R 2 . Furthermore, b(s, t) ∈ R 2 denotes a body force per unit reference length and the stretch ν(s, t) = ∂ s r is introduced as a measurement of longitudinal deformation of the string. The motion of the extensible string is governed by a quasilinear hyperbolic PDE (cf. [7]) which again fits with the square matrix B = E 2ρ 1 − 1 ν 2 I + 2 ν 4 (∂ s r ⊗ ∂ s r) and f (t) = F (t)ν(0, t)/N (ν(0, t)) into the framework of problem (1) by assuming the constitutive relation N (ν) = (EA/2ν)(ν 2 − 1). The numerical results (cf. Fig.2) have been achieved by applying bi-linear shape functions in space-time.