A variational formulation for motion design of adaptive compliant structures

Adaptive structures are characterized by their ability to adjust their geometrical and other properties to changing loads or requirements during service. This contribution deals with a method for the design of quasi-static motions of structures between two prescribed geometrical configurations that are optimal with regard to a specified quality function while taking large deformations into account. It is based on a variational formulation and the solution by two finite element discretizations, the spatial discretization (the standard finite element mesh) and an additional discretization of the deformation path or trajectory. For the investigations, an exemplary objective function, the minimization of the internal energy, integrated along the deformation path, is used. The method for motion design presented herein uses the Newton-Raphson method as a second order optimization algorithm and allows for analytical sensitivity analysis. The proposed method is verified and its properties are investigated by benchmark examples including rigid body motions, instability phenomena and determination of inextensible deformations of shells.


INTRODUCTION
Energy efficiency and sustainability play an increasing role in engineering and architecture. Reducing the amount of material used for construction does not only save resources but also reduces embedded energy. One possibility to realize extreme lightweight design is to make use of adaptive structures that optimally adjust their geometry to current and changing conditions by active motion. Here, two fundamentally different types of adaption via geometry change can be distinguished.
The first type of adaptive structure serves the purpose of adapting to changing loads or requirements and is mostly referred to as a smart structure. A system of sensors and actuators allows the structure to adjust, respond and actively counteract to varying loads. Thus, stresses, deformations or vibrations may be reduced in order to ensure usability and serviceability [1,2,3,4]. This approach potentially allows for eminent material savings. Examples for this type of adaptive structure are the Smart Shell [5] or an "infinitely stiff" cantilever beam [6] as well as a highrise building, which serves as a demonstrator in a current research project dealing with such adaptive structures [7]. The shape transition is realized by different types of actuators, from piezoelectric or electroactive polymer actuators [8,9,10], up to shape-memory alloys [11], which are usually controlled by optimal control algorithms [12]. The aforementioned types of adaptive structures are characterized by the fact that the individual states differ only by minor changes in geometry such that geometrically linear analyses are sufficient.
The second type of adaptive structure is not intended to adapt to varying loads but to changing requirements during service. This is, for example, the case for deployable and retractable structures. Prominent examples are the opening and closing of roofs, especially of stadiums (e.g. the Commerzbank-Arena in Frankfurt, Germany [13]) or adaptive folding bridges (e.g. the Kiel Hrn Footbridge in Kiel, Germany [14]). But also adaptive faade elements, beyond conventional sun-blinds, which can be opened and closed depending on the position of the sun and the control of interior daylight, belong to this kind of adaptive structures and may significantly contribute to the energy efficiency of a building [15]. Realisations are e.g. the One Ocean Expo 2012 Pavilion in Korea [16] or the biomimetic faade elements Flectofin [17] and Flectofold [18]. Another current research field, where this type of adaptive structure plays a prominent role, is the shape change of morphing wings of airplanes [19,20,21,22,23].
In this case, the geometries of the individual configurations differ significantly from each other. The standard approach to achieve variability in geometry is the targeted introduction of joints and hinges between stiff elements and associated defined kinematics by unfolding, sliding and similar mechanisms. But particularly these joints and hinges frequently represent weak spots of the structure and may be prone to failure. Another strategy for geometrical variability are discrete systems with integrated actuators. Such systems, which can also account for large deformations, can take the form of trusses [24,25,26], tensegrity structures [27,28,29,30,31] or lattice structures [32], where individual actuator elements can vary their length and therefore change the shape of the entire structure. This strategy is in contrast to an overall flexibility that is distributed in the entire structure and is therefore able to perform a smoothly distributed motion. This is, for example, the case for pure bending deformations, so-called inextensional deformations, in flexible and shape changing shells [33,34]. There are also approaches available to combine the discrete flexibility by joints with a distributed structural flexibility. Recent research investigates the optimization of flexibility and compliance of structures in order to enable an efficient deformation. These so-called compliant or morphing structures [35,36,37,38,39,40,41,42,43] are characterized by continuous stiffness changes and a varying stiffness distribution within the structure, leading to the formation of specific hinge zones. The challenge here is that, despite compliance, the structure remains strong enough to still be able to withstand loads in all configurations. The concept of multistable compliant structures also represents a possibility to deal with this problem and, at the same time, to keep the configurations stable without continuously expending effort [44,45].
Not only the geometries of the individual configurations at the beginning and the end of the motion have to meet specified requirements, but also the shape transition between the configurations. These transitions and movements imply stress onto the structure and require energy. The way of controlling the actuators plays a decisive role in the efficiency of the morphing structure.
In control theory, especially in optimal control and in motion planning of robots [46,47], exactly this problem of optimal trajectories is already addressed. Most structures that are investigated and calculated in these research areas, e.g. robots, are characterized by a discrete kinematic description and no (or negligible) elastic deformation. This leads to fewer degrees of freedom compared to continuously morphing structures. However, in the field of continuum robotics and hyper-redundant manipulation, exactly such continuous robots or systems are planned and investigated. A review of this field is given in Rus and Tolley [48]. But also in this case, consideration of a large number of degrees of freedom in combination with motion planning and optimal control strategies is challenging. Therefore, simplifications are made to capture the kinematics, like a piecewise constant curvature (PCC) model [49]. This enables control of such systems and solution of the inverse kinematics problem (calculating the required curvatures for a given end position), but planning and optimal control methods for continuous robots without a PCC model still remains a challenge and an open research tasks.
There are already examples and methods where the mechanics and analysis of structures are combined with optimal control and motion planning strategies, especially with regard to adaptive structures and large deformations. Ibrahimbegovic et al. [50] successfully combined optimal control with non-linear structural mechanics of a beam to reach a deformed end configuration with certain properties. Furthermore, Veuve et al. [30], Sychterz and Smith [31] and Masic and Skelton [29] used motion planning algorithms for the trajectory planning of tensegrity structures.
The focus of this study, however, is analysis and optimization with respect to certain requirements of the trajectory itself between prescribed configurations of any kind of flexible structure without a control algorithm and related sensoric equipment. In Werter et al. [51], the required actuation energy was already considered for the design of a morphing airplane wing and Maute and Reich [19] combined geometry optimization of a compliant wing structure with additional optimization of the adaption mechanism. Also, the pure transition between configurations can be formulated as an optimization problem. However, when non-linear kinematics is considered any evaluation of the objective function consists of a complete non-linear analysis, including all pertinent issues like convergence problems, instabilities and high computational cost. The use of higher order optimization methods requires fewer evaluations of the objective function but the required sensitivity analyses again cause computational expense.
This contribution deals with shape transitions as a motion between two (or more) geometrical configurations, taking geometrically non-linear structural behavior into account. The idea is to design shape transitions based on a variational formulation. A quasi-static process is assumed such that no inertia effects are considered. The solution is found by a finite element discretization of the path (trajectory) along with a Newton-Raphson solution algorithm, which can be interpreted as a second order optimization algorithm. Due to the path discretization, analytical sensitivities can be calculated by making use of standard components of the spatial finite elements, e.g. the stiffness matrix.
The paper is organized as follows. First, we refer to the Brachistochrone problem, which represents one of the first problems solved by variational principles, as an illustrative example to motivate the use of a variational method for motion design. Chapter 2 presents the problem statement and its solution with different strategies. In Section 3, a finite element solution algorithm is described to solve the problem of motion design of structures. The underlying weak form is based on a functional to be minimized to obtain an "optimal" path. The method is developed with the exemplary objective of minimizing the integral of the strain energy along the entire motion path. By a discretization of the motion path with finite elements, in addition to spatial discretization, a non-linear system of equations is obtained. Chapter 4 presents several numerical experiments to verify the proposed method by means of problems with known exact solutions, for instance, the motion of a kinematic system with zero strain energy throughout the entire process. Furthermore, motions dominated by instability behavior (following the motivation by [52,53]) and further potential of motion design and applications for a number of flexible structures, like the calculation of inextensible deformations of shells, are investigated. Finally, some conclusions are given and open issues and potential future developments are discussed in Chapter 5.

Historical background
To motivate the use of a variational formulation for the solution, the classical Brachistochrone problem is considered and similarities to the problem of motion design are highlighted. The Brachistochrone problem is one of the first problems that was solved with the calculus of variations and it represents the base for its development. In 1696, in the journal "Acta Eruditorum", published by Gottfried Wilhelm Leibniz, Johann Bernoulli posed to the scientific community the following challenge: Given two points A and B in a vertical plane, what is the curve traced out by a point acted on only by gravity, which starts at A and reaches B in the shortest time? [54] Shortly after the publication of the problem, Bernoulli received a letter from Leibniz in which he explained that he "is attracted by the problem like Eva by the apple", but at the same time, he asked for an extension of the deadline, since the problem reached other countries only after a few months. Bernoulli agreed with the proposal and reformulated the problem. At the beginning of 1697, an anonymous solution to the problem appeared in the journal "Philosophical Transactions of the Royal Society of London" [55] and finally in May 1697 Leibniz published a collection of the submitted solutions [56], in which also the anonymous solution, which Bernoulli identified directly as Newton's solution ("from the claw of the lion"), was reprinted. The solution by Jacob, Johann's brother, was then further developed and a few years later Leonhard Euler named it the calculus of variations [57].
where m represents the mass, v the velocity, y is the vertical abscissa and mgy A = const. defines the reference energy. Solving for the velocity yields v = 2g y A − y(x) . ( By using the definition of the velocity as time derivative of the arc length v = ds dt , an infinitesimal time increment can be written as dt = 1 v ds. The total time required for traveling from A to B is thus The infinitesimal arc length ds can be calculated from the Pythagorean theorem, An illustration for this derivation is given in Figure 1. Combining eq. (3) and eq. (4) yields the functional for the Brachistochrone problem

Euler Lagrange equation and exact solution
As the integrand F of the functional T (eq. 5) depends on the function y and its derivative y only and beyond this not explicitly on the variable x itself, the simplified Euler-Lagrange equation [58] F − y ∂F ∂y = C can be used. Application of eq. (6) to the functional of the Brachistochrone problem yields A substitution is used for solution of the problem, where trigonometric functions and a parametric form turn out to be a clever choice: The complete derivation is omitted at this point but can be looked up in the appendix in Chapter A. Eventually, a parametric representation of the coordinates x and yis obtained that is the function of a cycloid It must be noted thatt does neither represent the time nor the arc length, but is the angle by which a rolling circle has rotated, a point of which generates the curve (x(t), y(t)). The constants C 1 and C 2 as well as the parameter valuet E at point B are derived by the boundary conditions at the starting point A and the endpoint B: x(t = 0) = x A , x(t =t E ) = x B and y(t =t E ) = y B , which themselves represent non-linear functions that need to be solved iteratively. The condition y(t = 0) = y A is fulfilled by definition of the problem. An exemplary solution with fixed points A and B is given in Figure 2 (left).

Solution with finite elements
In general variational problems, a closed form solution is often not available. In such cases, an approximate solution may be obtained by the finite element method. This requires the functional to be formulated in a parametric form where s represents a path parameter. As the total length of the solution curve is initially unknown, the integration bound s B is unknown as well. Therefore, a mapping parameter su is defined, enabling integration over a variables and a fixed and specified domains ∈ [0, 1] The mapping parameter contains information about the arc length itself Inserting this into the functional yields Variation with respect to the unknown functions x and y then follows as This is the weak form of the Brachistochrone problem. The next steps follow the standard procedure of a finite element formulation. First, a discretization for x and y as well as their variations δx and δy is introduced where the matrix N contains the shape functions and the vectors x, y, δx and δy contain discrete nodal values of the unknowns x and y, respectively. In this case, the discretization is the same for every function, following the Bubnov-Galerkin approach. By inserting the discretization into the variation moving the vectors δx and δy out of the integral and applying the discrete form of the fundamental lemma of the calculus of variations, a residual is obtained. After linearization it can be solved iteratively for the nodal values x and y, which provide an approximation for the solution functions x and y in a parametric form. This discretization represents a discretization of the path that the point with mass m follows from A to B and is therefore referred to as path discretization in the following. Figure 2 (center) shows the result with a path discretization with linear Lagrange shape functions and 15 elements. As the parametrization of the curve, and thus the placement of the nodes along the solution curve, is not unique, an equal length of the path elements is enforced for the regularization of the solution. This extra constraint is enforced by Lagrange multipliers. It can be seen that the solution with finite elements and linear shape functions approximates well the exact curve, obtained in 2.2.2, but, due to the linear functions, it still contains kinks. The number of degrees of freedom is twice the number of internal nodes (x-and y-coordinate at each node). An improvement of the approximation is possible by using an interpolation by B-spline-functions, which can also be seen in Figure 2 (right), where the solution for a discretization with two elements and cubic shape functions is illustrated. The higher continuity enables a better approximation without kinks and fewer internal nodes (or control points) resulting in fewer degrees of freedom. This problem formulation of the Brachistochrone can be taken as a simplified template for what is intended to be done in motion design. The goal is finding a path between two configurations, e.g. the points A and B or an open and a closed geometry of an adaptive element that fulfills specified demands like minimization of the total required time, energy or effort for traversing from A to B, or any other objective.

MOTION DESIGN OF STRUCTURES
3.1. Functional and its variation as the starting point for motion design 3.1.1. Integrated internal energy as objective function The objective of motion design is mathematically expressed as minimization of a functional. In the Brachistochrone problem this was the required time for the mass point to run from A to B. In motion design, it can be any property that the motion is expected to have. One objective could be to minimize the straining that a structure is subjected to while undergoing a certain motion. This is comparable to the dimensionless quantity cost of transport that is used in various disciplines like biology and robotics. It represents a measure to quantify the cost or energy efficiency of different transport methods, i. e. walking, swimming or flying of an animal or driving of a vehicle from one location to an other. Here, the notion is transferred to a cost of deformation for flexible structures, where an energy criterion is utilized. Thus, the internal energy integrated over the entire motion path s is chosen here as an exemplary functional that represents the strain energy integrated along the motion path: Assuming small strains (but large displacements and rotations), a linear elastic St. Venant-Kirchhoff material law is used for the relationship between Green-Lagrange strain and second Piola-Kirchhoff stress. Incidentally, it is pointed out that this functional serves as a proof of concept and can be replaced by other objectives.
For further explanation of the problem and objective, an illustrating example is introduced in Figure 3. A simple truss structure (Young's modulus E = 30000, cross section area A = 0.1), forming a shallow arc, is supposed to deform from a starting configuration, shown in black, to a target configuration, shown in blue. This scenario is obviously inspired by the bi-stable setup of a snap-through problem. The blue target configuration, however, is not the stress-free snapped-through configuration of the black one but deviates from it by a horizontal shift of the central node.
A motion that minimizes the integrated internal energy, i.e. the functional, is to be found. The structure contains two unconstrained displacement degrees of freedom and it is assumed that forces can be applied to both of them to reach the target configuration. During deformation, the point P follows the, yet unknown, trajectory (red) until it arrives at the end position P . This trajectory can also be found on the plane which is spanned by the axes D 1 and D 2 in the diagram on the right in Figure 3. On the vertical axes, the internal energy Π int throughout the motion is plotted as a curved line, lying on the corresponding potential function. The area of the resulting surface is the value of the functional. The goal of this specific motion design task is to find the trajectory (the motion) that minimizes this area.

Specification of motion and arc length
The length of the path, along which the internal energy of the structure is integrated, can be associated with the arc length of the displacement field (cf. Figure 4) of the Figure 3. Two bar truss with initial configuration, target configuration and visualization of the functional underlying motion, which in turn is a function of the position X of the structure as well as the progress of the motion, the pseudo time t . (18) In order to consider the motion in its entirety within the functional, the internal energy is integrated along the deformation path s. This deformation path s represents a scalar measure that indicates by how much the structure has already moved and deformed. It is defined here as the arc length of the displacement field u(X, t). To obtain one scalar quantity from the displacement field, depending on the position vector X (Figure 4), the mean value of the displacement arc length inside the spatial domain Ω is used. Based on the same derivation as in the Brachistochrone problem (see Figure 1 and eq. (4), an infinitesimal arc length can then be specified for a three-dimensional problem as including the volume V of the domain. In the illustrating example of Figure 3, the arc length turns out to be the length of the trajectory of point P multiplied with the length of one bar (due to symmetry), because this problem involves only two degrees of freedom that are located at the same node. As this length of the trajectory is initially unknown, the integration limits of the functional are not fixed and remain unknown, as it was the case in the Brachistochrone problem. Therefore, another parameter must be introduced that indicates the motion progress with fixed integration bounds. In quasi-static structural analysis, often a normalized pseudo-time t is used, which runs from t = 0 to t = 1. This idea is adopted here and the motion parameter is re-defined as a normalized arc length of the deformation path. To use this path parameter as integration variable, a substitution is necessary, The mapping parameter is referred to as total arc length and the functional J transforms to 3.2. Spatial discretization and path discretization 3.2.1. General concept In order to solve the variational problem in equation (25), two discretizations are introduced. They divide both the spatial domain and the path into elements. Thus, a continuous problem is transferred into a discrete problem with a finite amount of degrees of freedom. Those degrees of freedom are located at the nodes forming the elements.
3.2.2. Spatial discretization First, a standard spatial discretization of Ω is introduced. The same continuity requirements apply as in a standard non-linear structural finite element analysis. The space is divided into n ele subdomains Ωe, the finite elements, on which integration is performed, The unknown displacement field is approximated by shape functions, interpolating the unknowns between discrete values at n nd,ele nodes per element With the mapping of Ωe to a reference element with the natural coordinates ξ by a Jacobian Je = ∂Xe ∂ξ the interpolation of the displacement field, as well as the reference and current geometry in an isoparametric concept, can be expressed as In a Bubnov-Galerkin approach, identical interpolation functions are used for the approximation of the variations In order to obtain a system of equations for all parameters of the problem, the nodal values of the single elements need to be assembled. This can be formally written with an assembly operator containing the information about element connectivity (topology). The nodal values of the displacement field and the current configuration still depend on the path variable s. The dimension of the vector D(s) is equal to the number of spatial degrees of freedom n dof .

Path discretization
As the spatial degrees of freedom D(s) are still functions of the path, a second discretization, the path discretization is required. It can also be denoted as a discretization of motion. This differs from a discretization in time, because the path also depends on the deformation of the structure, whereas time is considered an independent and autonomous value. The path, parametrized by the normalized arc lengths ∈ [0, 1], is subdivided inton ele path elements The shape functions can either be defined in the normalized parameter spaces ∈ [0, 1] or they can be transformed by a Jacobian. Variables referring to the path discretization are marked with a bar( •). Element numbers are indicated with a superscript (instead of a subscript, as in spatial discretization) for distinction. Also for the path elements, interpolation functions, a mapping to a reference element by a Jacobian Jē = ∂sē ∂ξ , as well as a Bubnov-Galerkin approach are used The nodes of the path discretizationk represent the different geometric configurations throughout the motion, including the initial, intermediate and end configurations, The shape functions serve for interpolation between the individual configurations and the total degrees of freedom are all n dof spatial degrees of freedom in every configurationk. Therefore, the vectord consists of where the length ofd is the amount of total degrees of freedom in one path element Path elements are one-dimensional. In the case of two spatial degrees of freedom, as in the illustrative example ( Figure 5), the path elements discretize the trajectory of point P . With an increasing number of degrees of freedom in space, they form a one-dimensional subspace within an n dof -dimensional hyperspace. The matrix of shape functions can contain any type of function. In Figure 5 a discretization with linear Lagrange shape functions is illustrated, but B-splines and higher order functions are also possible. As the variational index in the calculation of the arc length is equal to 1, at least C 0 -continuous functions are needed. Assembly is performed in the same manner as in spatial discretization with the assembly operator Path discretization is further visualized in the illustrative example of the two bar truss in Figure 5, where the vectors Dk are displayed. In this example, the vectors consist of two components due to the two spatial degrees of freedom For this prescribed path discretization by four linear elements, the degrees of freedom per path element arē The vector with all degrees of freedom can then be built by assemblȳ where the parameters D 0 are already defined in the problem formulation as the initial and the target configuration. The method for motion design therefore aims to move the intermediate configurations such that the functional J is minimized.
One issue with the motion design problem described so far is its potential ill-posedness for special cases. Within path discretization, nodes may be located anywhere on the trajectory, while still approximating the same curve (take a straight line as the simplest example for such a situation). Thus, the solution is no unique. This issue is well-known from, for instance, shape optimization and form finding problems of thinwalled structures, where nodes can be dislocated in-plane without changing the geometry. A corresponding regularization can be realized by either enforcing a constant path element size or by controlling the increments of a specified displacement degree of freedom throughout the deformation process. This aspect is further elaborated in Section 3.6.
At a first glance, path discretization resembles time integration in dynamic problems by space-time finite elements. However, both approaches are fundamentally different, since the arc length depends on the deformation of the structure, whereas time represents an independent and autonomous value. Another difference lies in the application of the two approaches. While space-time elements are mostly used to calculate and represent dynamic problems containing inertia effects, the motion path discretization is developed for quasi-static loading situations and static problems. This has an impact on the required element size, as dynamic effects, which can potentially be missed by using a time discretization that is too coarse, do not play a role in motion design problems.
In the next section, the discretized variations of the individual terms are introduced, where the two discretizations for space and path are introduced separately and successively.
The first variation of the Green-Lagrange strain with respect to the spatially discrete parameters d(s) can then be written as The strain-displacement matrix of the spatial elements is still continuous in the path variable s.
Total arc length The total arc length from eq. (23) is now expressed in a spatially discretized form. First, the lengths of the trajectories of the individual nodes of the spatial discretization are generated. The length of the nodal trajectories follows as A simple summation of the lengths of the nodal trajectories results in a dependency of the spatial discretization. This would mean that the total arc length of a motion of a coarse spatial discretization is smaller compared to the one of a finer mesh. Therefore, a mean value, in this case the root mean square, of the nodal trajectory lengths is determined. To calculate the average value, the influence volume V k of each individual node k is determined as illustrated in Figure 6. The root mean square can then be computed as This represents the spatially discretized total arc length. As the total arc length s s u depends on the derivative of the total displacements only, its first variation is which includes the gradient of su with respect to the derivatives of the displacement degrees of freedom. For a concise notation, the following abbreviations for the derivatives of the total arc length with respect to the spatial parameters are introduced As the vectors δd s and δd s ,s only contain discrete values in Ω, they can be extracted from the integral and by assembly eq. (55) can be written as where denotes the usual assembly operator. In this equation, the global vector of internal forces F s int and the internal energy Π s int , both still continuous in the path s, are identified Linearization Both terms can now be linearized separately Inserting these terms into the variation leads to the spatially discretized linearized variation This can further be modified by extraction of δD from the integral and by rearranging the terms to

Global linearized system of equations
From the condition that the discretized variation must vanish for any δd the following system of equations can be derived With the definitions we obtain the system of equations in the familiar format, Note that with this system the entire problem is solved monolithically, instead of incrementally proceeding along the path. On convergence of the iterative solution method, all intermediate configurations along the path are obtained in one go.
The system of equations depends on the used element as it includes the stiffness matrix and the internal forces. However, all ingredients can be combined in a modular manner. Thus, it does not pose any problem to use various element types, like mixed elements or isogeometric spatial discretizations.
In the solution of the illustrating two bar truss problem, the path is discretized by 14 linear path elements. The predictor represents an entire motion and it is chosen as a linear interpolation between the initial configuration and the target configuration.
The solution process of the non-linear problem with Newton's method converges after 9 iterations below the tolerance value of 10 −8 of the L 2 norm of the residual. As a result of this simple motion design problem, it is found, that for minimizing the integrated internal energy it is beneficial to first enforce a purely vertical snap-through, followed by a horizontal movement (as opposed to following the straight path of the linear Figure 7. Predictor motion and solution of illustrating example predictor motion). Therefore, the motion design method yields an optimized motion in a purely formalized way without the need to put any engineering expert knowledge into the analysis. The trajectory of the midpoint (as well as the path) is longer in the solution than in the predictor, but the proposed detour leads to a smaller accumulated internal energy throughout the motion. The difference is visualized in a plot of the internal energy over the two spatial degrees of freedom in Figure 7 (center) and a projection of the resulting functional surfaces where the internal energy is plotted versus the arc length in Figure 7 (right). The snapthrough characteristics can also be detected in the progress of the internal energy for the final solution. After snap-through, the internal energy reaches the value zero.
In order to realize the prescribed deformation that results from motion design, forces are needed. Those forces are evaluated after convergence from the internal forces and equilibrium of internal and external forces. This means that for the special case considered so far, where all degrees of freedom are controlled, the equilibrium conditions are not needed for the solution of the motion design problem, but the equilibrium equations can be used in post-processing of the nodal forces.

A generalized system of equations for any objective function
So far, the minimization of the internal energy, integrated along the path, was used as a proof of concept for the proposed motion design framework. However, in principle any functional or objective function can be used. In general, we define a quantity F that depends on the displacements and integrate it along the path to obtain the generalized functional After spatial discretization, the variation is built by the product rule The linearization can then be derived for the two terms Path discretization and reordering the terms leads tō (83) This is the system of equations for a general objective function for which analytical derivatives can be calculated. The minimized quantity must depend on the displacements to apply this method for motion design, but they can as well be calculated numerically. In a lot of cases, quantities and their derivatives, e.g. strains and stresses, can be used that are routinely available in standard finite element codes.

Aspects of convergence
The derived non-linear problem needs to be solved iteratively. The degrees of freedom are all spatial degrees of freedom in every single configuration. Some configurations are known, like the initial, starting geometry and the final, target geometry. It is also possible to define only parts of the target geometry. As the entire motion has to be found by one monolithic solution of the system of equations, the predictor describes an entire motion. It can be seen in Figures 3,5 and 7 that is a problem with two spatial degrees of freedom in one point P , the path elements discretize the trajectory of this point P until it reaches the endpoint P . However, the distribution and length of the path elements are not yet specified, which leads to an ill-posed problem. This can be fixed by additional controls. Either the progression of one (or multiple) spatial degrees of freedom can be prescribed by e.g. constant increments between the configurations or equal length of the path elements can be enforced by the introduction of Lagrange multipliers and the corresponding changes in the system of equations. The second method, the use of Lagrange multipliers, was applied in the illustrating example (see Figure 7). As the difference between the predictor motion and the final result may be large, the solution process sometimes suffers from convergence problems and the Newton method occasionally diverges after a number of iterations. There is no straightforward analogy to incremental-iterative solution procedures with the option to decrease the size of the increments in order to improve convergence behavior. There are, however, some other measures that can be taken: Fewer degrees of freedom by path approximation with B-splines The path may be approximated more efficiently by B-spline functions compared to linear Lagrange functions. The resulting reduction in the number of degrees of freedom leads to better convergence behavior.
Improved predictor from solution with coarse path discretization Likewise, a calculation with a small number of path elements, and therefore fewer degrees of freedom, improves the convergence of the Newton algorithm. This might result in a poor approximation due to the coarse discretization. This solution, however, can be used as an improved predictor for a computation with a finer path discretization. More generally, a hierarchically modified predictor improves convergence.
Better predictor by a preanalysis This method also focuses on the improvement of the first guess, the predictor. Instead of a linear interpolation, a standard non-linear analysis of the structure can be carried out to already approach a feasible motion. To this end, the internal forces in the end configurations are taken as the external forces (load case) for the non-linear analysis. The obtained equilibrium path is then used as a predictor motion.
Modification of the Newton method with a relaxation factor To improve convergence, a modification of the Newton method can be applied in which a relaxation factor prevents off-shooting from a possible solution during iterations in which the norm of the residual increases. This method has been presented in [59] and further investigated and developed in [60].
All the described methods can also be combined.

NUMERICAL EXPERIMENTS
Numerical experiments are presented to test and verify the potential of the proposed method for motion design. Some examples serve for quantitative benchmarking of the solution and others contribute to a better understanding of the solution and possible applications.

Kinematic structures for benchmarking
First, benchmarking examples, for which the exact solutions are known, are used for verification. Obvious scenarios are kinematic mechanisms for which the internal energy is identically zero throughout the entire motion.

Kinematic truss system
The first example is a kinematic truss system with four nodes, three bars and two supports, as shown in Figure 8. This kinematic system allows for a purely energy-free rigid body movement during which the lengths of the bars do not change. Pretending that the target configuration is unknown, it is sufficient to specify the vertical displacement of the second node to obtain a well-posed problem. The other displacements are expected to adjust to allow the kinematic movement (to minimize the functional of motion design). The path is discretized by 14 linear Lagrange elements, resulting inn dof = 42 degrees of freedom. For regularization, the vertical displacement of the second node is prescribed throughout the motion. The predictor is -intentionally naive -chosen to be a linear interpolation between the initial and the prescribed end configuration for the upper left node, while the upper right node does not move at all, as seen in Figure 8. As this is not a rigid body motion, forces are needed to enforce it, which are shown as red arrows. The predictor is far off the expected solution, with a functional value of J = 12843.
The table in Figure 8 shows a comparison of seven snapshots, i.e. every second intermediate configuration, of the converged solution and the predictor motion. It can be seen that the solution obtained from motion design provides the expected result with zero length changes of the individual bars despite the naive predictor. The difference between the linear interpolation and the final solution is also visible in the diagram on the right, where the internal energy is plotted versus the arc length along the path. The value of the functional represents the integral of the curve, i.e. the area of the blue and red area, respectively. The length Figure 8. Kinematic structure for benchmarking of the path differs between the two motions. In the predictor motion, the moving node moves directly to the end position while the other node does not move at all. This results in a shorter path compared to the resulting path of the solution. The internal energy is much smaller as almost no strain is present in the bars. The fact that the energy is not exactly zero results from the error due to path discretization errors (note the difference in the y-axis of the factor 10 4 ).
The value of J = 0.05 of the functional is not exactly zero for the obtained solution due to the error from path discretization with linear elements. By a refinement of the path discretization, as seen in Figure 9, center, the approximation quality increases and the value of the functional approaches zero. The analysis with a discretization by B-splines, shown on the right, enables an even better approximation of the curved motion trajectory and results in a smaller value of the functional with fewer degrees of freedom. Figure 9. Convergence study for kinematic structure 4.1.2. Folding motion with quadrilateral elements In the next example, a fold-like motion of an assembly of four quadrilateral elements, as shown in Figure 10, is modeled. The elements are connected with hinges either on the upper or on the lower corner. This enables mirroring of the geometry solely by rigid body translations and rotations. The path is approximated with quadratic B-splines andn ele = 6 elements. Again, the target geometry is assumed to be unknown, only the vertical displacement of the upper second and fourth node is prescribed and the vertical displacement increments are controlled during motion. The predictor motion is a linear interpolation and shows an unphysical movement with self-penetration of the elements.
By an iterative solution with the linearized system of equations presented in eq. (77), the correct motion with zero internal energy throughout the motion is found (see Figure 10).
The analysis of the two kinematic structures verifies that the proposed method finds the correct solution for this specific class of problems.

Motion design with multiple snap-through processes
The combination of three pairs of hinged bars, shown in Figure 11, represents a system for which the equilibrium path may exhibit multiple limit points, i.e. horizontal tangents, where snap-through occurs. The upper two-bar truss with a larger cross-sectional area A 2 is supported by two other two-bar trusses (cross section A 1 ) that can perform snap-through as well. The path is discretized by 32 linear elements. Only the vertical displacement of the upper node is prescribed and controlled throughout the motion. To enhance convergence behavior, the predictor is calculated and updated hierarchically from a solution obtained with a course path discretization as explained in Section 3.6. Therefore the linear interpolation with 32 elements does not represent the predictor motion in this case.
The solution is compared to a linear interpolation between the initial and a mirrored geometry, which is expected to represent a better approximation than the naive linear interpolation of only the upper central node to the target position. This results in a motion dominated by global snap-through. The result of motion design provides a different type of motion. When the side structures don't perform the snap-through at the same time, internal energy can be "saved" in the upper truss. This behavior can also be detected in the Figure 11. Motion design with a combination of multiple snap-through progress of the internal energy. The surface of the two "snap-through bulges" are clearly identifiable. The resulting end configuration is found to be the horizontally mirrored geometry, which reduces the internal energy back to zero. The value of the functional decreases significantly from J = 440 to J = 11.

Motion design in a bifurcation problem
In a high two-bar truss subject to a vertical load, as shown in Figure 12, bifurcation occurs before a limit point (snap-through) is reached, as it is the case in a shallow two-bar truss. Here, a system with a width-to-height ratio of 1 : 3 is investigated. The path is discretized by 20 linear elements and the vertical displacement is controlled for motion design. For the vertically flipped geometry as target configuration, linear interpolation describes a purely vertical snap-through motion. Indeed, this happens to represent a stationary point for the functional of motion design. However, it provides a relative maximum of J, not a minimum, meaning that it is a worst case scenario. Therefore, the predictor needs to be modified significantly to improve convergence of the motion design algorithm to the desired solution.
For example, instead of a linear interpolation, a combination of the primary path -up to the critical point -followed by an arbitrarily chosen branch of the secondary equilibrium path, describing the deformation after buckling of the structure, can be used as predictor. The optimized motion found on the basis of this predictor is shown in Figure 12 and yields a functional value of J = 779. It is significantly smaller than the value of J = 1449 obtained from linear interpolation, the "worst case scenario" mentioned above. But it is also superior to the value J = 845 obtained for the improved predictor based on the secondary path, which confirms the virtue of the method of motion design.
It can be observed, however, that the maximum value of the internal energy during deformation is higher for the optimized motion than for the secondary path (diagram on the right in Figure 12). The fact that the functional value is still lower for the optimized motion follows from two aspects: During the first phase of the deformation process, the internal energy value is higher in the predictor than in the optimized motion and the deformation path is slightly longer. These aspects are dominant and lead to the reduction of the functional value, even though the maximum value of internal energy is higher in the optimized solution.
Yet an alternative predictor is the so-called critical path. It is defined as the path that connects configurations for which the determinant of the stiffness matrix is zero, det K = 0. It leads to a functional value of J = 956, which is worse than both the optimal solution and the solution obtained from following the secondary path. Nevertheless, it is a valid predictor for obtaining convergence of the motion design algorithm. Figure 12. Analysis of a two-bar truss with bifurcation and motion design 4.2.4. Snap-through of a shallow arch The last example with a snap-through is a shallow arc, which is modeled as a two-dimensional structure under plane stress conditions, using quadrilateral finite elements, as shown in Figure 13. The fully prescribed target geometry is artificially chosen and represents the (approximately) mirrored geometry of the initial configuration. The path is discretized by five elements with cubic B-splines as shape functions and the vertical displacement of the center node is controlled for motion design.
First, purely displacement-based bilinear quadrilateral elements are used for spatial discretization. The predictor motion is again a linear interpolation between the initial and end configuration and represents a symmetric snap-through-dominated motion. By motion design, an antisymmetric swaying motion is found, which decreases the value of the functional from J = 1263 to J = 144. In this symmetric example, the mirrored deformation is equivalent to the calculated solution. It is well known that displacement-based finite elements suffer from locking. Therefore, the influence of locking on the result of motion design is investigated next by using quadrilateral finite elements including an Enhanced Assumed Strain (EAS) formulation, proposed by [61]. The resulting stiffness matrix and internal forces of this formulation can simply be plugged into the system of equations from Chapter 3. With four additional strain parameters per element, shear locking and volumetric locking can be eliminated. Already in the snapshots of the motion shown in Figure 13 the difference in the result obtained with finite elements that suffer from locking and locking-free elements is visible, although the overall character of the motion seems to be similar. Even though locking is not very dominant in this example, it can be observed that the EASelements exhibit more bending throughout the motion. The artificial energy that results from locking effects increases the internal energy along the path and acts as a penalty for bending modes. Locking-free elements avoid this penalty and the value of the functional decreases significantly from J = 144 to J = 58. It can therefore be expected that for structures that are more prone to locking, like slender thin-walled structures, locking has a significant effect on the result of motion design. Corresponding observations have been made for optimization problems in [62].

Specification of intermediate configurations
Beyond the possibility to specify initial and target configuration, also intermediate configurations can be included as an objective for motion design. Figure 14 shows a three-dimensional curved cantilever beam, discretized by trilinear volume elements. There are two intermediate and a final target configuration. First, the cantilever tip is rotated by -90 • around the z-axis (Configuration 1). Configuration 2 is defined as a straight, vertical bar. The final target configuration is identical to the initial configuration, rotated by 90 • about the z-axis, as shown in Figure 14, right.
Path discretization is accomplished with a total of nine quadratic elements -three elements for every deformation stage -using B-splines as shape functions. While path discretization with quadratic B-splines is usually C 1 -continuous, continuity is reduced to C 0 at path nodes that correspond to the intermediate configurations in order to respect the expected non-smoothness of the solution.
For stabilization, the displacement in y-direction of one node at the cantilever tip was controlled in stage 1. In this case, the C 1 -continuity of the path discretization was reduced to a C 0 -continuity at the node of configuration 1.
The final solution with the individual stages is shown in Figure 14.

Calculation of inextensible deformations of shells
4.4.1. Basic concept Motion design can also be performed for shells. One interesting option in this context is a modification of the functional by replacing the complete internal energy by the membrane energy only. This provides a method to compute motions that try to avoid membrane strains during deformation while bending remains without any penalization. The results are (nearly) inextensional deformations. Inextensional deformations of surfaces are defined as deformations that preserve lengths and angles of infinitesimal line elements at each point. Gaussian curvature remains constant during inextensional deformations. For thin shells (and beams) inextensional deformations can also be classified as pure bending deformations. In the following examples, isogeometric Kirchhoff-Love elements, as presented in [63], are used. It has to be noted that these elements still suffer from membrane locking, although by integration of the internal energy, strain oscillations are leveled out to a certain extent. However, when recovering the forces required to realize the found deformations, the effect of locking leads to values that are too large.

4.4.2.
Deformation of a cantilever beam One important special case are inextensional deformations of developable structures with Gaussian curvature equal to zero, e.g. bending of a cylinder to a flat plane. The deformation of a cantilever beam illustrated in Figure 15 represents the same phenomenon in a simple Figure 14. Specification of intermediate configurations on a cantilever beam modelled with volume elements two-dimensional configuration. The left side is clamped and for the target configuration, the final location of the tip is prescribed. It is defined in a way that allows the final configuration to be a perfect half-circle.
Initially, the beam is discretized with only two quadratic isogeometric elements to improve convergence due to the low number of degrees of freedom. The path is discretized by two quadratic elements with B-spline shape functions. By motion design, an inextensional deformation is found, where the straight cantilever is bent to a half-circle while preserving its length. Figure 15. Motion design of a cantilever with shell elements and corresponding inextensible deformations However, despite the good geometry approximation by using NURBS as shape functions, this mesh is too coarse to provide reasonable results in terms of stress and strain (and therefore the internal energy). Therefore, in the following an improved approximation of the geometry is realized by using 12 quadratic elements and the motion obtained with the coarse mesh is used as a predictor. The resulting motion is shown in Figure 15 and it closely resembles the one obtained with the coarse mesh, but does not represent the perfect half circle. The value of the functional, however, significantly decreases from J = 4.8 to J = 0.01.
It has to be noted that in this numerical experiment the solution is not unique. Any deformed geometry for which the cantilever has the same length as the original flat configuration can be reached by an inextensional deformation. Accordingly, the problem is ill-posed. Nevertheless, one valid solution is found with apparently no numerical problems. In order to understand this surprising phenomenon one has to first understand that the finite elements used are based on a displacement based standard Galerkin formulation with no measures to avoid locking. For the problem at hand, membrane locking is crucial. For the given discretization with 12 quadratic elements the effect is not very strong. However, the corresponding parasitic non-zero membrane strains are large enough to have a regularizing effect on the process of motion design.

4.4.3.
Transformation of a helicoid to a catenoid A classical example for an inextensional deformation is the transformation of a helicoid to a catenoid, shown in Figure 16. It is a rare example from the field of analytical differential geometry for which an analytical solution for large inextensional deformations exists in the case of Gaussian curvature being non-zero.
The helicoid is discretized with 4 × 4 cubic elements with B-spline shape functions. For the target geometry, only the final position of the upper and lower Edge (A − B,C − D) are prescribed. The path is also coarsely discretized with two quadratic elements with B-spline shape functions. For motion design, the vertical displacement of a point at the upper edge is controlled.
Since the final geometry is only defined at the edges, the predictor, again obtained from a linear interpolation, shows a relatively bad first guess. The solution of the motion design problem not only determines the correct inextensible deformation, but also the correct final geometry, the catenoid. The value of the functional is J = 0.2, i.e. again close to zero. In this paper, a variational method for the design of motions of continuously deformable structures between two geometrical configurations, fulfilling certain desired properties and considering large displacements, has been presented. As a proof of concept, a functional defining the integrated internal energy along the motion path was defined. A combination of spatial discretization and path discretization is used for the numerical solution of the underlying problem. The convergence behavior of the motion design problem is enhanced by various methods. In problems, where continuous deformation paths are expected, B-spline shape functions can be used to reduce the number of degrees of freedom compared to standard Lagrange discretization. An advantageous side effect of fewer degrees of freedom is a further improvement of convergence behavior in the iterative process.
The applied solution procedure can be interpreted as a second order optimization algorithm. For the problems studied herein, the derivatives (sensitivities) are calculated analytically. In the given framework this can be accomplished for any kind of finite element for spatial discretization.
Implementation of the motion design method was successfully verified by benchmark problems for which the exact solution is known. Additionally, the effect of instability and snap-through phenomena for motions with the prescribed functional was investigated in corresponding examples. The feasibility of the method for the design of inextensible deformations of shells was also demonstrated.
The successful application of motion design to detect or develop kinematic mechanisms reveals a genuine potential for application to adaptive and deployable structures. The evolution of the required actuation forces during the deformation is recovered after the motion is found. The restriction to those cases results from the fact, that the realization of the designed motion potentially requires forces at every degree of freedom, whereas usually, prescribed load cases or actuators are at hand. The future and subsequent steps in the development of the motion design method are therefore investigations on how to incorporate discrete actuator elements and that a resulting motion can be realized only with a restricted and specified amount of possible load cases.

A. EXACT SOLUTION OF THE BRACHISTOCHRONE PROBLEM
The starting point is eq.
Solving it for y yields y = 1 2gC 2 1 (y A − y) A substitution is necessary for a solution and a clever choice for this problem is the parametric representation of trigonometrical functions