A variational formulation for motion design of adaptive compliant structures
AA variational formulation for motion design of adaptive compliantstructures
Renate Sachse , ∗ and Manfred Bischoff Institute for Structural Mechanics, University of Stuttgart,Pfaffenwaldring 7, D-70550 Stuttgart, Germany
SUMMARYAdaptive 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 designof quasi-static motions of structures between two prescribed geometrical configurations that are optimalwith regard to a specified quality function while taking large deformations into account. It is based on avariational formulation and the solution by two finite element discretizations, the spatial discretization (thestandard finite element mesh) and an additional discretization of the deformation path or trajectory. For theinvestigations, an exemplary objective function, the minimization of the internal energy, integrated alongthe deformation path, is used. The method for motion design presented herein uses the Newton-Raphsonmethod as a second order optimization algorithm and allows for analytical sensitivity analysis. The proposedmethod is verified and its properties are investigated by benchmark examples including rigid body motions,instability phenomena and determination of inextensible deformations of shells.KEY WORDS: Optimization, adaptive structures, deployable structures, motion design, variationalformulation, shape morphing structure, compliant mechanism1. INTRODUCTIONEnergy efficiency and sustainability play an increasing role in engineering and architecture. Reducing theamount 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 optimallyadjust their geometry to current and changing conditions by active motion. Here, two fundamentally differenttypes of adaption via geometry change can be distinguished.The first type of adaptive structure serves the purpose of adapting to changing loads or requirements andis 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 reducedin order to ensure usability and serviceability [1, 2, 3, 4]. This approach potentially allows for eminentmaterial 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 projectdealing 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], whichare usually controlled by optimal control algorithms [12]. The aforementioned types of adaptive structuresare characterized by the fact that the individual states differ only by minor changes in geometry such thatgeometrically linear analyses are sufficient.The second type of adaptive structure is not intended to adapt to varying loads but to changingrequirements 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 ∗ Correspondence to: Renate Sachse, Institute for Structural Mechanics, University of Stuttgart, Pfaffenwaldring 7, D-70550 Stuttgart, Germany, E-mail: [email protected] a r X i v : . [ c s . C E ] J u l R. SACHSE AND M. BISCHOFF and closed depending on the position of the sun and the control of interior daylight, belong to this kind ofadaptive structures and may significantly contribute to the energy efficiency of a building [15]. Realisationsare 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 prominentrole, 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. Thestandard approach to achieve variability in geometry is the targeted introduction of joints and hingesbetween 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 proneto failure. Another strategy for geometrical variability are discrete systems with integrated actuators. Suchsystems, which can also account for large deformations, can take the form of trusses [24, 25, 26], tensegritystructures [27, 28, 29, 30, 31] or lattice structures [32], where individual actuator elements can vary theirlength and therefore change the shape of the entire structure. This strategy is in contrast to an overallflexibility that is distributed in the entire structure and is therefore able to perform a smoothly distributedmotion. 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 discreteflexibility by joints with a distributed structural flexibility. Recent research investigates the optimization offlexibility and compliance of structures in order to enable an efficient deformation. These so-called compliantor 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. Thechallenge here is that, despite compliance, the structure remains strong enough to still be able to withstandloads in all configurations. The concept of multistable compliant structures also represents a possibilityto deal with this problem and, at the same time, to keep the configurations stable without continuouslyexpending effort [44, 45].Not only the geometries of the individual configurations at the beginning and the end of the motion haveto meet specified requirements, but also the shape transition between the configurations. These transitionsand movements imply stress onto the structure and require energy. The way of controlling the actuators playsa 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 thisproblem of optimal trajectories is already addressed. Most structures that are investigated and calculated inthese 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 continuousrobots or systems are planned and investigated. A review of this field is given in Rus and Tolley [48]. Butalso in this case, consideration of a large number of degrees of freedom in combination with motion planningand 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 ofthe inverse kinematics problem (calculating the required curvatures for a given end position), but planningand optimal control methods for continuous robots without a PCC model still remains a challenge and anopen research tasks.There are already examples and methods where the mechanics and analysis of structures are combinedwith optimal control and motion planning strategies, especially with regard to adaptive structures and largedeformations. Ibrahimbegovic et al. [50] successfully combined optimal control with non-linear structuralmechanics of a beam to reach a deformed end configuration with certain properties. Furthermore, Veuveet al. [30], Sychterz and Smith [31] and Masic and Skelton [29] used motion planning algorithms for thetrajectory planning of tensegrity structures.The focus of this study, however, is analysis and optimization with respect to certain requirements ofthe trajectory itself between prescribed configurations of any kind of flexible structure without a controlalgorithm and related sensoric equipment. In Werter et al. [51], the required actuation energy was alreadyconsidered for the design of a morphing airplane wing and Maute and Reich [19] combined geometryoptimization 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, whennon-linear kinematics is considered any evaluation of the objective function consists of a complete non-linearanalysis, 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 therequired sensitivity analyses again cause computational expense.
This contribution deals with shape transitions as a motion between two (or more) geometricalconfigurations, taking geometrically non-linear structural behavior into account. The idea is to design shapetransitions based on a variational formulation. A quasi-static process is assumed such that no inertia effectsare considered. The solution is found by a finite element discretization of the path (trajectory) along with aNewton-Raphson solution algorithm, which can be interpreted as a second order optimization algorithm. Due
VARIATIONAL FORMULATION FOR MOTION DESIGN OF STRUCTURES
2. THE BRACHISTOCHRONE PROBLEM
To motivate the use of a variational formulation for the solution, the classical Brachistochrone problem isconsidered and similarities to the problem of motion design are highlighted. The Brachistochrone problemis one of the first problems that was solved with the calculus of variations and it represents the base for itsdevelopment. In 1696, in the journal “Acta Eruditorum”, published by Gottfried Wilhelm Leibniz, JohannBernoulli 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 ononly 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 explainedthat he “is attracted by the problem like Eva by the apple”, but at the same time, he asked for an extension ofthe deadline, since the problem reached other countries only after a few months. Bernoulli agreed with theproposal and reformulated the problem. At the beginning of 1697, an anonymous solution to the problemappeared in the journal “Philosophical Transactions of the Royal Society of London” [55] and finally in May1697 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. Thesolution by Jacob, Johann’s brother, was then further developed and a few years later Leonhard Euler namedit the calculus of variations [57].
Starting point for the calculation is conservation ofenergy with the kinetic energy E kin and the potential energy E pot E kin + E pot = 12 mv + mgy ( x ) = mgy A , (1)where m represents the mass, v the velocity, y is the vertical abscissa and mgy A = const . defines thereference energy. Solving for the velocity yields v = (cid:113) g (cid:0) y A − y ( x ) (cid:1) . (2)By using the definition of the velocity as time derivative of the arc length v = d s d t , an infinitesimal time increment can be written as d t = v d s . The total time required for traveling from A to B is thus T = (cid:90) T d t = (cid:90) s B s A (cid:113) g (cid:0) y A − y ( x ) (cid:1) d s = min . (3) R. SACHSE AND M. BISCHOFF
The infinitesimal arc length d s can be calculated from the Pythagorean theorem, d s = (cid:112) d x + d y = (cid:112) d x + d y d x d x = (cid:114)(cid:16) d x d x (cid:17) + (cid:16) d y d x (cid:17) d x = (cid:112) y (cid:48) ( x ) d x. (4)An illustration for this derivation is given in Figure 1. Combining eq. (3) and eq. (4) yields the functionalfor the Brachistochrone problem T = (cid:90) s B s A (cid:113) g (cid:0) y A − y ( x ) (cid:1) d s = (cid:90) x B x A (cid:115) y (cid:48) ( x ) g (cid:0) y A − y ( x ) (cid:1) d x. (5) Figure 1. Brachistochrone problem
As the integrand F of the functional T (eq. 5) dependson the function y and its derivative y (cid:48) only and beyond this not explicitly on the variable x itself, thesimplified Euler-Lagrange equation [58] F − y (cid:48) ∂F∂y (cid:48) = C (6)can be used. Application of eq. (6) to the functional of the Brachistochrone problem yields (cid:115) y (cid:48) g ( y A − y ) − y (cid:48) y (cid:48) (cid:112) g ( y A − y )(1 + y (cid:48) ) = 1 (cid:112) g ( y A − y )(1 + y (cid:48) ) = C . (7)A substitution is used for solution of the problem, where trigonometric functions and a parametric form turnout to be a clever choice: gC ( y A − y ) = sin (cid:18) ¯ t (cid:19) = 12 (cid:0) − cos(¯ t ) (cid:1) . (8)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 y is obtained that is the function of a cycloid y (¯ t ) = y A − gC (cid:0) − cos(¯ t ) (cid:1) x (¯ t ) = 14 gC (¯ t − sin(¯ t )) + C . (9)It must be noted that ¯ t does neither represent the time nor the arc length, but is the angle by which a rollingcircle has rotated, a point of which generates the curve ( x (¯ t ) , y (¯ t )) . The constants C and C as well asthe parameter value ¯ t E at point B are derived by the boundary conditions at the starting point A and theendpoint B: x (¯ t = 0) = x A , x (¯ t = ¯ t E ) = x B and y (¯ t = ¯ t E ) = y B , which themselves represent non-linearfunctions that need to be solved iteratively. The condition y (¯ t = 0) = y A is fulfilled by definition of theproblem. An exemplary solution with fixed points A and B is given in Figure 2 (left). 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 requiresthe functional to be formulated in a parametric form T = (cid:90) s B s A =0 (cid:113) g (cid:0) y A − y ( s ) (cid:1) d s, (10) VARIATIONAL FORMULATION FOR MOTION DESIGN OF STRUCTURES s represents a path parameter. As the total length of the solution curve is initially unknown, theintegration bound s B is unknown as well. Therefore, a mapping parameter s u is defined, enabling integrationover a variable ¯ s and a fixed and specified domain ¯ s ∈ [0 , (cid:90) s B ( . . . ) d s = (cid:90) ( . . . ) d s d¯ s d¯ s = (cid:90) ( . . . ) s u d¯ s. (11)The mapping parameter contains information about the arc length itself s u := d s d¯ s = (cid:112) d x + d y d¯ s = (cid:115)(cid:18) d x d¯ s (cid:19) + (cid:18) d y d¯ s (cid:19) = (cid:112) x (cid:48) (¯ s ) + y (cid:48) (¯ s ) . (12)Inserting this into the functional yields T = (cid:90) (cid:115) x (cid:48) (¯ s ) + y (cid:48) (¯ s ) g (cid:0) y A − y (¯ s ) (cid:1) d¯ s (13)Variation with respect to the unknown functions x and y then follows as δT = (cid:90) (cid:32)(cid:115) x (cid:48) + y (cid:48) g (cid:0) y A − y (cid:1) δy + x (cid:48) δx (cid:48) + y (cid:48) δy (cid:48) (cid:113) g (cid:0) y A − y (cid:1)(cid:0) x (cid:48) + y (cid:48) (cid:1) δx (cid:48) (cid:33) d¯ s = 0 . (14)This is the weak form of the Brachistochrone problem. The next steps follow the standard procedure of afinite element formulation. First, a discretization for x and y as well as their variations δx and δy is introduced x ≈ x h = Nx δx ≈ δx h = N δ x y ≈ y h = Ny δy ≈ δy h = N δ y (15)where the matrix N contains the shape functions and the vectors x , y , δ x and δ y contain discrete nodalvalues 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 δT = (cid:90) (cid:32)(cid:115) ( N (cid:48) x ) + ( N (cid:48) y ) g (cid:0) y A − ( Ny ) (cid:1) N δ y + N (cid:48) xN (cid:48) δ x + N (cid:48) yN (cid:48) δ y (cid:113) g (cid:0) y A − ( Ny )) (cid:1)(cid:0) ( N (cid:48) x ) + ( N (cid:48) y ) (cid:1) (cid:33) d¯ s = 0 , (16)moving the vectors δ x and δ y out of the integral and applying the discrete form of the fundamental lemmaof the calculus of variations, a residual is obtained. After linearization it can be solved iteratively for thenodal values x and y , which provide an approximation for the solution functions x and y in a parametricform. This discretization represents a discretization of the path that the point with mass m follows from Ato 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 15elements. 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. Thisextra constraint is enforced by Lagrange multipliers. It can be seen that the solution with finite elements andlinear 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 Figure 2. Exact solution (left), solution with finite elements by a linear (center) and a cubic B-splinediscretization (right)
R. SACHSE AND M. BISCHOFF
B-spline-functions, which can also be seen in Figure 2 (right), where the solution for a discretization withtwo elements and cubic shape functions is illustrated. The higher continuity enables a better approximationwithout 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 isintended to be done in motion design. The goal is finding a path between two configurations, e.g. the pointsA and B or an open and a closed geometry of an adaptive element that fulfills specified demands likeminimization of the total required time, energy or effort for traversing from A to B, or any other objective.3. MOTION DESIGN OF STRUCTURES
The objective of motion design is mathematicallyexpressed as minimization of a functional. In the Brachistochrone problem this was the required time forthe mass point to run from A to B. In motion design, it can be any property that the motion is expected tohave. One objective could be to minimize the straining that a structure is subjected to while undergoing acertain motion. This is comparable to the dimensionless quantity cost of transport that is used in variousdisciplines 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 onelocation to an other. Here, the notion is transferred to a cost of deformation for flexible structures, where anenergy criterion is utilized. Thus, the internal energy integrated over the entire motion path s is chosen hereas an exemplary functional that represents the strain energy integrated along the motion path: J = (cid:90) s Π int d s = (cid:90) s (cid:90) Ω E T S dΩ d s = (cid:90) s (cid:90) Ω E T CE dΩ d s = min . (17)Assuming small strains (but large displacements and rotations), a linear elastic St. Venant-Kirchhoffmaterial 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 otherobjectives.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 . ), forming a shallow arc,is supposed to deform from a starting configuration, shown in black, to a target configuration, shown inblue. This scenario is obviously inspired by the bi-stable setup of a snap-through problem. The blue targetconfiguration, however, is not the stress-free snapped-through configuration of the black one but deviatesfrom 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 structurecontains two unconstrained displacement degrees of freedom and it is assumed that forces can be applied toboth 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 (cid:48) . This trajectory can also be found on the plane whichis spanned by the axes D and D in the diagram on the right in Figure 3. On the vertical axes, the internalenergy Π 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 taskis to find the trajectory (the motion) that minimizes this area. The length of the path, along which the internal energy of thestructure 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
VARIATIONAL FORMULATION FOR MOTION DESIGN OF STRUCTURES X of the structure as well as the progress ofthe motion, the pseudo time t u ( X , t ) = (cid:34) u ( X , t ) u ( X , t ) u ( X , t ) (cid:35) . (18)In order to consider the motion in its entirety within the functional, the internal energy is integrated alongthe deformation path s . This deformation path s represents a scalar measure that indicates by how muchthe 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 thesame derivation as in the Brachistochrone problem (see Figure 1 and eq. (4), an infinitesimal arc length canthen be specified for a three-dimensional problem as d s = 1 V (cid:90) Ω (cid:113) d u + d u + d u dΩ (19) including the volume V of the domain. In the illustrating example of Figure 3, the arc length turns out tobe the length of the trajectory of point P multiplied with the length of one bar (due to symmetry), becausethis problem involves only two degrees of freedom that are located at the same node. As this length of thetrajectory is initially unknown, the integration limits of the functional are not fixed and remain unknown, asit was the case in the Brachistochrone problem.Therefore, another parameter must be introduced that indicates the motion progress with fixed integrationbounds. 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 thedeformation path. To use this path parameter as integration variable, a substitution is necessary, (cid:90) s u ( . . . ) d s = (cid:90) ( . . . ) d s d¯ s d¯ s = (cid:90) ( . . . ) s u d¯ s. (20)The mapping parameter s u := d s d¯ s = V (cid:82) Ω (cid:112) d u + d u + d u dΩd¯ s (21) = 1 V (cid:90) Ω (cid:115)(cid:18) d u d¯ s (cid:19) + (cid:18) d u d¯ s (cid:19) + (cid:18) d u d¯ s (cid:19) dΩ (22) = 1 V (cid:90) Ω (cid:113) u , ¯ s + u , ¯ s + u , ¯ s dΩ (23)is referred to as total arc length and the functional J transforms to J = (cid:90) (cid:90) Ω E T CE dΩ s u d¯ s = min . (24) Figure 4. Illustration of the displacement field
R. SACHSE AND M. BISCHOFF
The strain and the total arc length are functions of the unknowndisplacements. The variation is computed according to the chain rule and set equal to zero, δJ = (cid:90) (cid:20) (cid:90) Ω δ E T CE dΩ s u + (cid:90) Ω E T CE dΩ δs u (cid:21) d¯ s = 0 . (25) In order to solve the variational problem in equation (25), two discretizations areintroduced. They divide both the spatial domain and the path into elements. Thus, a continuous problem istransferred into a discrete problem with a finite amount of degrees of freedom. Those degrees of freedomare located at the nodes forming the elements.
First, a standard spatial discretization of Ω is introduced. The same continuityrequirements 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, (cid:90) Ω ( . . . ) dΩ = n ele (cid:88) e =1 (cid:90) Ω e ( . . . ) dΩ e . (26) The unknown displacement field is approximated by shape functions, interpolating the unknowns betweendiscrete values at n nd , ele nodes per element u ( X , s ) ≈ u h ( X , s ) = n nd , ele (cid:88) k =1 N k ( X ) d k ( s ) = N ( X ) d ( s ) . (27)With the mapping of Ω e to a reference element with the natural coordinates ξ by a Jacobian J e = ∂ X e ∂ ξ theinterpolation of the displacement field, as well as the reference and current geometry in an isoparametricconcept, can be expressed as u h ( ξ , s ) = n nd , ele (cid:88) k =1 N k ( ξ ) d k ( s ) = N ( ξ ) d ( s ) , (28) X h ( ξ , s ) = n nd , ele (cid:88) k =1 N k ( ξ ) X k = N ( ξ ) X , (29) x h ( ξ , s ) = n nd , ele (cid:88) k =1 N k ( ξ ) x k ( s ) = N ( ξ ) x ( s ) . (30)In a Bubnov-Galerkin approach, identical interpolation functions are used for the approximation of thevariations δ u ( ξ , s ) ≈ δ u h ( ξ , s ) = n nd , ele (cid:88) k =1 N k ( ξ ) δ d k ( s ) = N ( ξ ) δ d ( s ) . (31)In order to obtain a system of equations for all parameters of the problem, the nodal values of the singleelements need to be assembled. This can be formally written with an assembly operator D ( s ) = n ele (cid:91) e =1 d e ( s ) , (32)containing the information about element connectivity (topology). The nodal values of the displacementfield and the current configuration still depend on the path variable s . The dimension of the vector D ( s ) isequal to the number of spatial degrees of freedom n dof . VARIATIONAL FORMULATION FOR MOTION DESIGN OF STRUCTURES As the spatial degrees of freedom D ( s ) are still functions of the path, a seconddiscretization, the path discretization is required. It can also be denoted as a discretization of motion. Thisdiffers 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 normalizedarc length ¯ s ∈ [0 , , is subdivided into ¯ n ele path elements (cid:90) ( . . . ) d¯ s = ¯ n ele (cid:88) ¯ e =1 (cid:90) s ¯ e ( . . . ) d¯ s ¯ e . (33)The shape functions can either be defined in the normalized parameter space ¯ s ∈ [0 , or they can betransformed by a Jacobian. Variables referring to the path discretization are marked with a bar ¯( ◦ ) . Elementnumbers 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 ¯ e = ∂ ¯ s ¯ e ∂ ¯ ξ , as well as a Bubnov-Galerkin approach are used d ¯ h (¯ ξ ) = ¯ n nd , ele (cid:88) ¯ k =1 ¯ N ¯ k (¯ ξ )¯ d ¯ k = ¯ N ( ξ )¯ d , (34) δ d ¯ h (¯ ξ ) = ¯ n nd , ele (cid:88) ¯ k =1 ¯ N ¯ k (¯ ξ ) δ ¯ d ¯ k = ¯ N ( ξ ) δ ¯ d . (35)The nodes of the path discretization ¯ k represent the different geometric configurations throughout themotion, including the initial, intermediate and end configurations, ¯ d ¯ k = D ¯ k = D ( s = s ¯ k ) (36) δ ¯ d ¯ k = δ D ¯ k = δ D ( s = s ¯ k ) . (37)The shape functions serve for interpolation between the individual configurations and the total degrees offreedom are all n dof spatial degrees of freedom in every configuration ¯ k . Therefore, the vector ¯ d consists of ¯ n nd , ele subvectors ¯ d = (cid:104) D D . . . D ¯ k . . . D ¯ n nd , ele (cid:105) T (38) δ ¯ d = (cid:104) δ D δ D . . . δ D ¯ k . . . δ D ¯ n nd , ele (cid:105) T , (39)where the length of ¯ d is the amount of total degrees of freedom in one path element ¯ n dof = ¯ n nd , ele · n dof . (40)Path elements are one-dimensional. In the case of two spatial degrees of freedom, as in the illustrativeexample (Figure 5), the path elements discretize the trajectory of point P . With an increasing number ofdegrees of freedom in space, they form a one-dimensional subspace within an n dof -dimensional hyperspace. Figure 5. Illustrating example with motion elements and the discretized illustration of the functional0
R. SACHSE AND M. BISCHOFF
The matrix of shape functions can contain any type of function. In Figure 5 a discretization with linearLagrange shape functions is illustrated, but B-splines and higher order functions are also possible. As thevariational index in the calculation of the arc length is equal to 1, at least C -continuous functions areneeded. Assembly is performed in the same manner as in spatial discretization with the assembly operator ¯ D = ¯ n ele (cid:91) ¯ e =1 ¯ d ¯ e . (41)Path discretization is further visualized in the illustrative example of the two bar truss in Figure 5, wherethe vectors D ¯ k are displayed. In this example, the vectors consist of two components due to the two spatialdegrees of freedom D ¯ k = (cid:104) D ¯ k D ¯ k (cid:105) T . (42)For this prescribed path discretization by four linear elements, the degrees of freedom per path element are ¯ d = (cid:2) D D D D (cid:3) T (43) ¯ d = (cid:2) D D D D (cid:3) T (44) ¯ d = (cid:2) D D D D (cid:3) T (45) ¯ d = (cid:2) D D D end1 D end2 (cid:3) T . (46)The vector with all degrees of freedom can then be built by assembly ¯ D = (cid:2) D D D D D D D D D end1 D end2 (cid:3) T , (47)where the parameters D = 0 , D = 0 , D end1 and D end2 are already defined in the problem formulation as theinitial and the target configuration. The method for motion design therefore aims to move the intermediateconfigurations 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 thesame curve (take a straight line as the simplest example for such a situation). Thus, the solution is nounique. This issue is well-known from, for instance, shape optimization and form finding problems of thin-walled structures, where nodes can be dislocated in-plane without changing the geometry. A correspondingregularization can be realized by either enforcing a constant path element size or by controlling theincrements of a specified displacement degree of freedom throughout the deformation process. This aspectis further elaborated in Section 3.6.At a first glance, path discretization resembles time integration in dynamic problems by space-timefinite elements. However, both approaches are fundamentally different, since the arc length depends onthe deformation of the structure, whereas time represents an independent and autonomous value. Anotherdifference lies in the application of the two approaches. While space-time elements are mostly used tocalculate and represent dynamic problems containing inertia effects, the motion path discretization isdeveloped for quasi-static loading situations and static problems. This has an impact on the required elementsize, as dynamic effects, which can potentially be missed by using a time discretization that is too coarse, donot play a role in motion design problems.In the next section, the discretized variations of the individual terms are introduced, where the twodiscretizations for space and path are introduced separately and successively. Green-Lagrange strains
For improved readability, in the following parameters depending on the pathvariable s are written with a superscript s (cid:16) • (cid:17) ( s ) =: (cid:16) • (cid:17) s . (48) VARIATIONAL FORMULATION FOR MOTION DESIGN OF STRUCTURES d ( s ) canthen be written as δ E s = (cid:18) ∂ E s ∂ d s (cid:19) T δ d s . (49)The strain-displacement matrix B s = (cid:18) ∂ E s ∂ d s (cid:19) T (50)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 lengthof the nodal trajectories follows as s s u ,k = (cid:118) (cid:117) (cid:117)(cid:116) n disp,nd (cid:88) i D i,k ( s ) . (51)A simple summation of the lengths of the nodal trajectories results in a dependency of the spatialdiscretization. This would mean that the total arc length of a motion of a coarse spatial discretization issmaller compared to the one of a finer mesh. Therefore, a mean value, in this case the root mean square, ofthe nodal trajectory lengths is determined. To calculate the average value, the influence volume V k of eachindividual node k is determined as illustrated in Figure 6. The root mean square can then be computed as s s u = (cid:118)(cid:117)(cid:117)(cid:116) V n nd (cid:88) k V k s ,k . (52)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 δs s u = (cid:18) ∂s s u ∂ D s,s (cid:19) T δ D s,s , (53)which includes the gradient of s u 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 respectto the spatial parameters are introduced s s u := ∂s s u ∂ D s,s , S s u := ∂ s s u ( ∂ D s,s ) . (54) Figure 6. Illustration of the spatially discretized total arc length2
R. SACHSE AND M. BISCHOFF
Variation
Together with eq. (25) the discretized variation of the functional reads δJ = (cid:90) (cid:20) n ele (cid:88) e =1 (cid:90) Ω e ( δ d s ) T (cid:18) ∂ E s ∂ d s (cid:19) CE s dΩ e s s u + ( δ d s,s ) T s s u n ele (cid:88) e =1 (cid:90) Ω e E s T CE s dΩ e (cid:21) d¯ s = 0 . (55)As the vectors δ d s and δ d s,s only contain discrete values in Ω , they can be extracted from the integral andby assembly eq. (55) can be written as δJ = (cid:90) (cid:20) ( δ D s ) T n ele (cid:91) e =1 (cid:90) Ω e ( B s ) T S s dΩ e (cid:124) (cid:123)(cid:122) (cid:125) F s int s s u + ( δ D s,s ) T s s u n ele (cid:91) e =1 (cid:90) Ω e E s T S s dΩ e (cid:124) (cid:123)(cid:122) (cid:125) Π s int (cid:21) d¯ s = 0 , (56)where (cid:83) denotes the usual assembly operator. In this equation, the global vector of internal forces F s int andthe internal energy Π s int , both still continuous in the path s , are identified δJ = (cid:90) (cid:20) ( δ D s ) T F s int s s u + ( δ D s,s ) T s s u Π s int (cid:21) d¯ s = 0 (57) Linearization
Both terms can now be linearized separately
LIN ( F s int s s u ) = F s int s s u + ∂ F s int s s u ∂ D s ∆ D s + ∂ F s int s s u ∂ D s,s ∆ D s,s (58) = F s int s s u + (cid:18) ∂ F s int ∂ D s (cid:124) (cid:123)(cid:122) (cid:125) = K s T s s u + F s int ∂s s u ∂ D s (cid:124) (cid:123)(cid:122) (cid:125) =0 (cid:19) ∆ D s + (cid:18) ∂ F s int ∂ D s,s (cid:124) (cid:123)(cid:122) (cid:125) =0 s s u + F s int ∂s s u ∂ D s,s (cid:124) (cid:123)(cid:122) (cid:125) s s u (cid:19) ∆ D s,s (59) = F s int s s u + K s T s s u ∆ D s + F s int s s u ∆ D s,s (60)(61) LIN ( s s u Π s int ) = s s u Π s int + ∂ s s u Π s int ∂ D s ∆ D s + ∂ s s u Π s int ∂ D s,s ∆ D s,s (62) = s s u Π s int + (cid:18) ∂ s s u ∂ D s (cid:124) (cid:123)(cid:122) (cid:125) =0 Π s int + s s u ∂ Π s int ∂ D s (cid:124) (cid:123)(cid:122) (cid:125) = F s int (cid:19) ∆ D s + (cid:18) ∂ s s u ∂ D s,s (cid:124) (cid:123)(cid:122) (cid:125) S s u Π s int + s s u ∂ Π s int ∂ D s,s (cid:124) (cid:123)(cid:122) (cid:125) =0 (cid:19) ∆ D s,s ) (63) = s s u Π s int + s s u F s int ∆ D s + S s u Π s int ∆ D s,s . (64)Inserting these terms into the variation leads to the spatially discretized linearized variation δJ = (cid:90) (cid:20) ( δ D s ) T (cid:16) F s int s s u + K s T s s u ∆ D s + F s int s s u ∆ D s,s (cid:17) (65) + ( δ D s,s ) T (cid:16) s s u Π s int + s s u F s int ∆ D s + S s u Π s int ∆ D s,s (cid:17)(cid:21) d¯ s = 0 . (66) Eq. (66) is still continuous along the path, so the path discretization fromSection 3.2.3 is introduced for the parameters d s , their variation δ d s , the linearized parameters ∆ d s and thecorresponding partial derivatives D s = ¯ N ¯ D δ D s = ¯ N δ ¯ D ∆ D s = ¯ N ∆ ¯ D (67) D s,s = ¯ N ,s ¯ D δ D s,s = ¯ N ,s δ ¯ D ∆ D s,s = ¯ N ,s ∆ ¯ D . (68)Inserting the path discretization into eq. (66) yields the completely discretized and linearized variation δJ = ¯ n ele (cid:88) ¯ e =1 (cid:90) s ¯ e (cid:20) δ ¯ D T ¯ N T (cid:16) F int s u + K T s u ¯ N ∆ ¯ D + F int s u ¯ N ,s ∆ ¯ D (cid:17) (69) + δ ¯ D T ¯ N T ,s (cid:16) s u Π int + s u F int ¯ N ∆ ¯ D + S u Π int ¯ N ,s ∆ ¯ D (cid:17)(cid:21) d¯ s ¯ e = 0 . (70) VARIATIONAL FORMULATION FOR MOTION DESIGN OF STRUCTURES δ ¯ D from the integral and by rearranging the terms to δJ = ¯ n ele (cid:88) ¯ e =1 δ ¯ D T (cid:20) (cid:90) ¯ s ¯ e (cid:16) ¯ N T F int s u + ¯ N T ,s s u Π int (cid:17) d s ¯ e + (71) (cid:90) s ¯ e (cid:16) ¯ N T K T s u ¯ N + ¯ N T F int s u ¯ N ,s + ¯ N T ,s s u F int ¯ N + ¯ N T ,s S u Π int ¯ N ,s (cid:1) d s ∆ ¯ D (cid:21) d¯ s ¯ e = 0 (72) From the condition that the discretized variation must vanish for any δ ¯ d the following system of equationscan be derived ¯ n ele (cid:91) ¯ e =1 (cid:90) ¯ s ¯ e (cid:16) ¯ N T K T s u ¯ N + ¯ N T F int s u ¯ N ,s + ¯ N T ,s s u F int ¯ N + ¯ N T ,s S u Π int ¯ N ,s (cid:17) d s ¯ e ∆ ¯ D (73) = − ¯ n ele (cid:91) ¯ e =1 (cid:90) s ¯ e (cid:16) ¯ N T F int s u + ¯ N T ,s s u Π int (cid:17) d¯ s ¯ e . (74)With the definitions K mod = ¯ n ele (cid:91) ¯ e =1 (cid:90) ¯ s ¯ e (cid:16) ¯ N T K T s u ¯ N + ¯ N T F int s u ¯ N ,s + ¯ N T ,s s u F int ¯ N + ¯ N T ,s S u Π int ¯ N ,s (cid:17) d¯ s ¯ e (75) R mod = ¯ n ele (cid:91) ¯ e =1 (cid:90) ¯ s ¯ e (cid:16) ¯ N T F int s u + ¯ N T ,s s u Π int (cid:17) d¯ s ¯ e (76)we obtain the system of equations in the familiar format, K mod ∆ ¯ D = − R mod . (77)Note that with this system the entire problem is solved monolithically, instead of incrementally proceedingalong the path. On convergence of the iterative solution method, all intermediate configurations along thepath are obtained in one go.The system of equations depends on the used element as it includes the stiffness matrix and the internalforces. However, all ingredients can be combined in a modular manner. Thus, it does not pose any problemto 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 initialconfiguration and the target configuration.The solution process of the non-linear problem with Newton’s method converges after 9 iterations belowthe tolerance value of − of the L 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 verticalsnap-through, followed by a horizontal movement (as opposed to following the straight path of the linear Figure 7. Predictor motion and solution of illustrating example4
R. SACHSE AND M. BISCHOFF predictor motion). Therefore, the motion design method yields an optimized motion in a purely formalizedway without the need to put any engineering expert knowledge into the analysis. The trajectory of themidpoint (as well as the path) is longer in the solution than in the predictor, but the proposed detour leadsto a smaller accumulated internal energy throughout the motion. The difference is visualized in a plot of theinternal energy over the two spatial degrees of freedom in Figure 7 (center) and a projection of the resultingfunctional surfaces where the internal energy is plotted versus the arc length in Figure 7 (right). The snap-through characteristics can also be detected in the progress of the internal energy for the final solution. Aftersnap-through, the internal energy reaches the value zero.In order to realize the prescribed deformation that results from motion design, forces are needed. Thoseforces are evaluated after convergence from the internal forces and equilibrium of internal and externalforces. 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 equilibriumequations can be used in post-processing of the nodal forces.
So far, the minimization of the internal energy, integrated along the path, was used as a proof of concept forthe proposed motion design framework. However, in principle any functional or objective function can beused. In general, we define a quantity F that depends on the displacements and integrate it along the path to obtain the generalized functional J = (cid:90) F s u d¯ s = min . (78)After spatial discretization, the variation is built by the product rule δJ = (cid:90) (cid:104) δ ( D s ) T F s, D s s u + ( δ D s,s ) T s s u F s (cid:105) d¯ s = 0 . (79)The linearization can then be derived for the two terms LIN (cid:0) F s, D s s u (cid:1) = F s, D s s u + (cid:0) F s, D s D s s s u (cid:1) ∆ D s + (cid:0) F s, D s s s u (cid:1) ∆ D s,s (80) LIN (cid:0) s s u F s (cid:1) = (cid:0) s s u F s (cid:1) + (cid:0) s s u F s, D s (cid:1) ∆ D s + (cid:0) S s u F s (cid:1) ∆ D s,s . (81)Path discretization and reordering the terms leads to ¯ n ele (cid:88) ¯ e =1 (cid:90) ¯ s ¯ e (cid:16) ¯ N F , D s D s s u ¯ N + ¯ N F , D s s u ¯ N ,s + ¯ N T ,s s u F , D s ¯ N + ¯ N T ,s S u F ¯ N ,s (cid:17) d s ¯ e ∆ ¯ D (82) = − ¯ n ele (cid:88) ¯ e =1 (cid:90) s ¯ e (cid:16) ¯ N T F , D s s u + ¯ N T ,s s u F (cid:17) d¯ s ¯ e . (83)This is the system of equations for a general objective function for which analytical derivatives can becalculated. The minimized quantity must depend on the displacements to apply this method for motiondesign, 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. The derived non-linear problem needs to be solved iteratively. The degrees of freedom are all spatial degreesof freedom in every single configuration. Some configurations are known, like the initial, starting geometryand the final, target geometry. It is also possible to define only parts of the target geometry. As the entiremotion 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 freedomin one point P , the path elements discretize the trajectory of this point P until it reaches the endpoint P (cid:48) .However, the distribution and length of the path elements are not yet specified, which leads to an ill-posedproblem. This can be fixed by additional controls. Either the progression of one (or multiple) spatial degreesof freedom can be prescribed by e.g. constant increments between the configurations or equal length of thepath elements can be enforced by the introduction of Lagrange multipliers and the corresponding changes in VARIATIONAL FORMULATION FOR MOTION DESIGN OF STRUCTURES
Fewer degrees of freedom by path approximation with B-splines
The path may be approximated moreefficiently by B-spline functions compared to linear Lagrange functions. The resulting reduction inthe number of degrees of freedom leads to better convergence behavior.
Improved predictor from solution with coarse path discretization
Likewise, a calculation with a smallnumber of path elements, and therefore fewer degrees of freedom, improves the convergence ofthe 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 pathdiscretization. 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 carriedout to already approach a feasible motion. To this end, the internal forces in the end configurationsare taken as the external forces (load case) for the non-linear analysis. The obtained equilibrium pathis then used as a predictor motion.
Modification of the Newton method with a relaxation factor
To improve convergence, a modification ofthe Newton method can be applied in which a relaxation factor prevents off-shooting from a possiblesolution during iterations in which the norm of the residual increases. This method has been presentedin [59] and further investigated and developed in [60].All the described methods can also be combined.4. NUMERICAL EXPERIMENTSNumerical experiments are presented to test and verify the potential of the proposed method for motiondesign. Some examples serve for quantitative benchmarking of the solution and others contribute to a betterunderstanding of the solution and possible applications.
First, benchmarking examples, for which the exact solutions are known, are used for verification. Obviousscenarios are kinematic mechanisms for which the internal energy is identically zero throughout the entiremotion.
The first example is a kinematic truss system with four nodes, three barsand two supports, as shown in Figure 8. This kinematic system allows for a purely energy-free rigid bodymovement during which the lengths of the bars do not change. Pretending that the target configuration isunknown, it is sufficient to specify the vertical displacement of the second node to obtain a well-posedproblem. The other displacements are expected to adjust to allow the kinematic movement (to minimize thefunctional of motion design).The path is discretized by 14 linear Lagrange elements, resulting in ¯ n dof = 42 degrees of freedom.For regularization, the vertical displacement of the second node is prescribed throughout the motion. Thepredictor is – intentionally naive – chosen to be a linear interpolation between the initial and the prescribedend 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. Thepredictor 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 motiondesign provides the expected result with zero length changes of the individual bars despite the naivepredictor. The difference between the linear interpolation and the final solution is also visible in the diagramon the right, where the internal energy is plotted versus the arc length along the path. The value of thefunctional represents the integral of the curve, i.e. the area of the blue and red area, respectively. The length6 R. SACHSE AND M. BISCHOFF
Figure 8. Kinematic structure for benchmarking of the path differs between the two motions. In the predictor motion, the moving node moves directly tothe end position while the other node does not move at all. This results in a shorter path compared to theresulting 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 thedifference in the y-axis of the factor ).The value of J = 0 . of the functional is not exactly zero for the obtained solution due to the errorfrom 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 analysiswith a discretization by B-splines, shown on the right, enables an even better approximation of the curvedmotion trajectory and results in a smaller value of the functional with fewer degrees of freedom.Figure 9. Convergence study for kinematic structure In the next example, a fold-like motion of an assemblyof four quadrilateral elements, as shown in Figure 10, is modeled. The elements are connected with hingeseither on the upper or on the lower corner. This enables mirroring of the geometry solely by rigid bodytranslations and rotations. The path is approximated with quadratic B-splines and ¯ n 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 predictormotion 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 motionwith 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 solutionfor this specific class of problems. VARIATIONAL FORMULATION FOR MOTION DESIGN OF STRUCTURES
Another interesting aspect to further understand and validate the properties of theproposed motion design method is the analysis of structures and motions where snap-through or bifurcationcan occur. Various structures with potential instabilities are investigated next. A two-bar truss system thatperforms snap-through has already been presented during the derivation of the method in Chapter 3.
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-sectionalarea A is supported by two other two-bar trusses (cross section A ) that can perform snap-through as well.The path is discretized by 32 linear elements. Only the vertical displacement of the upper node is prescribedand controlled throughout the motion. To enhance convergence behavior, the predictor is calculated andupdated 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 isexpected to represent a better approximation than the naive linear interpolation of only the upper centralnode to the target position. This results in a motion dominated by global snap-through. The result of motiondesign provides a different type of motion. When the side structures don’t perform the snap-through at thesame 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-through8
R. SACHSE AND M. BISCHOFF progress of the internal energy. The surface of the two “snap-through bulges” are clearly identifiable. Theresulting end configuration is found to be the horizontally mirrored geometry, which reduces the internalenergy back to zero. The value of the functional decreases significantly from J = 440 to J = 11 . In a high two-bar truss subject to a vertical load, as shownin Figure 12, bifurcation occurs before a limit point (snap-through) is reached, as it is the case in a shallowtwo-bar truss. Here, a system with a width-to-height ratio of is investigated. The path is discretized by20 linear elements and the vertical displacement is controlled for motion design. For the vertically flippedgeometry 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 providesa relative maximum of J , not a minimum, meaning that it is a worst case scenario. Therefore, the predictorneeds to be modified significantly to improve convergence of the motion design algorithm to the desiredsolution.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 deformationafter buckling of the structure, can be used as predictor. The optimized motion found on the basis of thispredictor is shown in Figure 12 and yields a functional value of J = 779 . It is significantly smaller than thevalue 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, whichconfirms the virtue of the method of motion design.It can be observed, however, that the maximum value of the internal energy during deformation is higherfor the optimized motion than for the secondary path (diagram on the right in Figure 12). The fact that thefunctional value is still lower for the optimized motion follows from two aspects: During the first phase ofthe deformation process, the internal energy value is higher in the predictor than in the optimized motionand the deformation path is slightly longer. These aspects are dominant and lead to the reduction of thefunctional 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 connectsconfigurations for which the determinant of the stiffness matrix is zero, det K = 0 . It leads to a functionalvalue of J = 956 , which is worse than both the optimal solution and the solution obtained from followingthe secondary path. Nevertheless, it is a valid predictor for obtaining convergence of the motion designalgorithm. Figure 12. Analysis of a two-bar truss with bifurcation and motion design The last example with a snap-through is a shallow arc, which ismodeled 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 elementswith 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. Thepredictor motion is again a linear interpolation between the initial and end configuration and represents asymmetric 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, themirrored deformation is equivalent to the calculated solution. VARIATIONAL FORMULATION FOR MOTION DESIGN OF STRUCTURES to be similar. Even though locking is not very dominant in this example, it can be observed that the EAS-elements exhibit more bending throughout the motion. The artificial energy that results from locking effectsincreases the internal energy along the path and acts as a penalty for bending modes. Locking-free elementsavoid this penalty and the value of the functional decreases significantly from J = 144 to J = 58 . It cantherefore 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 madefor optimization problems in [62]. Beyond the possibility to specify initial and target configuration, also intermediate configurations can beincluded 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 astraight, 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 everydeformation stage – using B-splines as shape functions. While path discretization with quadratic B-splinesis usually C -continuous, continuity is reduced to C at path nodes that correspond to the intermediateconfigurations 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 -continuity of the path discretization was reduced to a C -continuity at the node ofconfiguration 1.The final solution with the individual stages is shown in Figure 14. Motion design can also be performed for shells. One interesting option in this contextis 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 whilebending remains without any penalization. The results are (nearly) inextensional deformations. Inextensionaldeformations of surfaces are defined as deformations that preserve lengths and angles of infinitesimal lineelements at each point. Gaussian curvature remains constant during inextensional deformations. For thinshells (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 hasto 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 requiredto realize the found deformations, the effect of locking leads to values that are too large.
One important special case are inextensional deformations ofdevelopable 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 simple0
R. SACHSE AND M. BISCHOFF
Figure 14. Specification of intermediate configurations on a cantilever beam modelled with volume elementstwo-dimensional configuration. The left side is clamped and for the target configuration, the final locationof 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 convergencedue to the low number of degrees of freedom. The path is discretized by two quadratic elements with B-splineshape functions. By motion design, an inextensional deformation is found, where the straight cantilever isbent to a half-circle while preserving its length.
Figure 15. Motion design of a cantilever with shell elements and corresponding inextensible deformations
VARIATIONAL FORMULATION FOR MOTION DESIGN OF STRUCTURES J = 4 . to J = 0 . .It has to be noted that in this numerical experiment the solution is not unique. Any deformed geometry forwhich the cantilever has the same length as the original flat configuration can be reached by an inextensionaldeformation. Accordingly, the problem is ill-posed. Nevertheless, one valid solution is found with apparentlyno numerical problems. In order to understand this surprising phenomenon one has to first understand thatthe finite elements used are based on a displacement based standard Galerkin formulation with no measuresto avoid locking. For the problem at hand, membrane locking is crucial. For the given discretization with 12quadratic elements the effect is not very strong. However, the corresponding parasitic non-zero membranestrains are large enough to have a regularizing effect on the process of motion design. A classical example for an inextensional deformation isthe transformation of a helicoid to a catenoid, shown in Figure 16. It is a rare example from the field ofanalytical differential geometry for which an analytical solution for large inextensional deformations existsin the case of Gaussian curvature being non-zero.The helicoid is discretized with × 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 pathis 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 linearinterpolation, shows a relatively bad first guess. The solution of the motion design problem not onlydetermines the correct inextensible deformation, but also the correct final geometry, the catenoid. The valueof the functional is J = 0 . , i.e. again close to zero.Figure 16. Transformation from a helicoid to a catenoid with a motion design analysis5. CONCLUSIONSIn this paper, a variational method for the design of motions of continuously deformable structures betweentwo geometrical configurations, fulfilling certain desired properties and considering large displacements, hasbeen 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 numericalsolution of the underlying problem. The convergence behavior of the motion design problem is enhanced byvarious methods. In problems, where continuous deformation paths are expected, B-spline shape functionscan be used to reduce the number of degrees of freedom compared to standard Lagrange discretization. Anadvantageous side effect of fewer degrees of freedom is a further improvement of convergence behavior inthe iterative process.2 R. SACHSE AND M. BISCHOFF
The applied solution procedure can be interpreted as a second order optimization algorithm. For theproblems studied herein, the derivatives (sensitivities) are calculated analytically. In the given frameworkthis 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 whichthe exact solution is known. Additionally, the effect of instability and snap-through phenomena for motionswith the prescribed functional was investigated in corresponding examples. The feasibility of the method forthe design of inextensible deformations of shells was also demonstrated.The successful application of motion design to detect or develop kinematic mechanisms reveals a genuinepotential for application to adaptive and deployable structures. The evolution of the required actuationforces during the deformation is recovered after the motion is found. The restriction to those cases resultsfrom the fact, that the realization of the designed motion potentially requires forces at every degree offreedom, whereas usually, prescribed load cases or actuators are at hand. The future and subsequent stepsin the development of the motion design method are therefore investigations on how to incorporate discreteactuator elements and that a resulting motion can be realized only with a restricted and specified amount ofpossible load cases. ACKNOWLEDGEMENTS
This work has been funded by the German Research Foundation (DFG) as part of the TransregionalCollaborative Research Centre (SFB/Transregio) 141 Biological Design and Integrative Structures/projectA04 under grant number INST 41/910-1 and the collaborative project Bio-inspirierte Materialsystemeand Verbundkomponenten fr nachhaltiges Bauen im 21ten Jahrhundert (BioElast), which is part of theZukunftsoffensive IV Innovation und Exzellenz Aufbau und Strkung der Forschungsinfrastruktur imBereich der Mikro- und Nanotechnologie sowie der neuen Materialien, funded by the State Ministry ofBaden-Wuerttemberg for Sciences, Research and Arts.
REFERENCES1. Housner GW, Bergman LA, Caughey TK, et al. Structural Control: Past, Present, and Future.
Journal of EngineeringMechanics
Journal of Structural Engineering
Computers & Structures
Journal of the International Association for Shell and Spatial Structures
Smart Materials and Structures
Steel Construction
EngineeringStructures
Sensors (Basel, Switzerland)
Mechanics of Soft Materials
Materials & Design (1980-2015)
Vibration Control of Active Structures: An Introduction . Solid Mechanics and Its ApplicationsSpringerNetherlands. 3 ed. 2011.13. G¨oppert K, Stein M. A Spoked Wheel Structure for the Worlds largest Convertible Roof The New CommerzbankArena in Frankfurt, Germany.
Structural Engineering International
Structural EngineeringInternational
Smart Materials and Structures
Bautechnik
Bioinspiration & Biomimetics
18. K¨orner A, Born L, Mader A, et al. Flectofolda biomimetic compliant shading device for complex free form facades.
Smart Materials and Structures
Journal of Aircraft
WIT Transactions on State ofthe Art in Science and Engineering . 2. WIT Press. 1 ed. 2006 (pp. 400–419)21. Vasista S, Tong L, Wong KC. Realization of Morphing Wings: A Multidisciplinary Challenge.
Journal of Aircraft
Journal of Aircraft
AerospaceScience and Technology
Robotics and Automation in Construction
InTech. 200825. Sofla AYN, Elzey DM, Wadley HNG. Shape morphing hinged truss structures.
Smart Materials and Structures
Structural andMultidisciplinary Optimization
Robotics andAutonomous Systems
Journal ofGuidance, Control, and Dynamics
Engineering Structures
Frontiers in Built Environment
YBL Journal of Built Environment
Bioinspiration & Biomimetics
Extremely Deformable Structures .562. Vienna: Springer Vienna. 2015 (pp. 179–267)35. Sigmund O. On the Design of Compliant Mechanisms Using Topology Optimization.
Mechanics of Structures andMachines
AIAA Journal
AnalogIntegrated Circuits and Signal Processing
Journal of Intelligent MaterialSystems and Structures
Journal of Intelligent Material Systems and Structures
Journal of Intelligent MaterialSystems and Structures
Smart Materials andStructures
Finite Elements in Analysis and Design
Proceedings of the 7th GACM Colloquium on ComputationalMechanics for Young Scientists from Academia and Industry
Institute for Structural Mechanics, University ofStuttgart; 2017: 61–65.44. Oh YS, Kota S. Synthesis of Multistable Equilibrium Compliant Mechanisms Using Combinations of BistableMechanisms.
Journal of Mechanical Design
International Journal of Solids and Structures
IEEE Access
Planning Algorithms . Cambridge: Cambridge University Press . 200648. Rus D, Tolley MT. Design, fabrication and control of soft robots.
Nature
TheInternational Journal of Robotics Research
International Journal for Numerical Methods in Engineering R. SACHSE AND M. BISCHOFF52. Reis PM. A Perspective on the Revival of Structural (In)Stability With Novel Opportunities for Function: FromBuckliphobia to Buckliphilia.
Journal of Applied Mechanics
Smart Materials and Structures
Acta Eruditorum
PhilosophicalTransactions of the Royal Society of London
Acta Eruditorum
Variationsrechnung . BI-Hochschultaschenb¨ucherMannheim: Bibliographisches Institut Mannheim .1970.59. Albanese R, Rubinacci G. Numerical procedures for the solution of nonlinear electromagnetic problems.
IEEETransactions on Magnetics
IEEE Transactions on Magnetics
InternationalJournal for Numerical Methods in Engineering
Computers & Structures
Computer Methods in Applied Mechanics and Engineering
A. EXACT SOLUTION OF THE BRACHISTOCHRONE PROBLEMThe starting point is eq. (5) in the main text, where the simplified Euler-Lagrange equation is applied for thefunctional of the Brachistochrone problem (cid:115) y (cid:48) g ( y A − y ) − y (cid:48) y (cid:48) (cid:112) g ( y A − y )(1 + y (cid:48) ) = 1 (cid:112) g ( y A − y )(1 + y (cid:48) ) = C . (84)Solving it for y (cid:48) yields y (cid:48) = (cid:115) gC ( y A − y ) − (cid:115) − gC ( y A − y )2 gC ( y A − y ) . (85)A substitution is necessary for a solution and a clever choice for this problem is the parametric representationof trigonometrical functions gC ( y A − y ) = sin (cid:18) ¯ t (cid:19) = 12 (cid:0) − cos(¯ t ) (cid:1) (86)and the resulting equation for yy = y A − gC sin (cid:18) ¯ t (cid:19) = y A − gC (cid:0) − cos(¯ t ) (cid:1) (87)is now the ansatz. Inserting it into eq. (85) y (cid:48) = (cid:118) (cid:117)(cid:117)(cid:117)(cid:116) − sin (cid:16) ¯ t (cid:17) sin (cid:16) ¯ t (cid:17) = cos (cid:16) ¯ t (cid:17) sin (cid:16) ¯ t (cid:17) = d y d x (88) VARIATIONAL FORMULATION FOR MOTION DESIGN OF STRUCTURES d x = sin (cid:16) ¯ t (cid:17) cos (cid:16) ¯ t (cid:17) d y. (89)The derivative of y with respect to ¯ t is still necessary d y d¯ t = 12 gC sin (cid:18) ¯ t (cid:19) cos (cid:18) ¯ t (cid:19) → d y = 12 gC sin (cid:18) ¯ t (cid:19) cos (cid:18) ¯ t (cid:19) d¯ t (90)and can now be inserted into eq. (89) d x = sin (cid:16) ¯ t (cid:17) cos (cid:16) ¯ t (cid:17) d y = sin (cid:16) ¯ t (cid:17) cos (cid:16) ¯ t (cid:17) gC sin (cid:18) ¯ t (cid:19) cos (cid:18) ¯ t (cid:19) d t = 12 gC sin (cid:18) ¯ t (cid:19) d t = 14 gC (cid:0) − cos(¯ t ) (cid:1) d¯ t. (91) By integration x = (cid:90) d x + C = (cid:90) gC (cid:0) − cos(¯ t ) (cid:1) d¯ t + C = 14 gC (¯ t − sin(¯ t )) + C (92)the equations for x and y can be obtained y = y A − gC sin (cid:18) ¯ t (cid:19) = y A − gC (cid:0) − cos(¯ t ) (cid:1) (93) x = 14 gC (¯ t − sin(¯ t )) + C . (94)As already explained in the main text, the constants C , C as well as the parameter value ¯ t E at pointB can be derived by the boundary conditions of the starting point A and the end point B: x ( t = 0) = x A , x (¯ t = ¯ t E ) = x B and y (¯ t = ¯ t E ) = y B , which represent themselves non-linear functions that need to be solvediteratively. The condition y (¯ t = 0) = y A is fulfilled by definition of the problem. Exemplary values aregiven for the starting point x A = 1 . , y A = 5 . and the end point x B = 10 . , y B = 2 . (calculated with therounded gravitation constant g = 10 ) C = 0 . C = 1 . t E = 4 .05