K C Park
email: kcpark@colorado.edu

Carlos A Felippa

Roger Ohayon

The d'Alembert-Lagrange principal equations and applications to floating flexible systems

Keywords: dynamics, d'Alembert's principle, floating systems, principal equations, Lagrange multipliers, partitioned analysis

This paper addresses the dynamics and quasi-statics of floating flexible structures as well as extensions to unconstrained substructures and partitions of coupled mechanical systems. The principal solution is defined as the state of self-equilibrated forces obtained as the particular solution of the rigid motion and interface equilibrium equations. This solution is independent of the stress-strain constitutive properties as well as of the compatibility equations. For statically determinate systems, the principal solution is the final force solution. For statically indeterminate systems, the correction due to flexibility and compatibility is orthogonal to the principal solution. The formulation is done in the context of d'Alembert's principle, which supplies the d'Alembert-Lagrange principal equations for floating bodies. These are obtained by summation of virtually working forces and moments acting on the floating systems. Applications of this approach are demonstrated on a set of dynamic and quasi-static example problems of increasing generality. Linkage to variational principles with an interface potential is eventually discussed as providing the theoretical foundation for handling interacting semi-discrete subsystems linked by node-collocated Lagrange multipliers.

INTRODUCTION

In presentations dealing with motions of fully unconstrained fluid and solid systems, the authors have often been asked whether a floating, flexible body under dynamic loads is in an equilibrium state, or whether an unsupported structure such as a flying aircraft maintains equilibrium as it moves and deforms. Doubts on such a fundamental subject suggest a gap in the teaching of classical The treatment of multiplier-connected components of a mechanical system has been addressed by many authors, particularly in the context of multibody dynamics [START_REF] Shabana | Dynamics of Multibody Systems[END_REF], parallel computations [START_REF] Farhat | A method of finite-element tearing and interconnecting and its parallel solution algorithm[END_REF][START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF], and multiphysics simulations [START_REF] Park | Partitioned formulation of internal fluid-structure interaction problems via localized Lagrange multipliers[END_REF][START_REF] Ross | Coupling and simulation of acoustic fluid-structure interaction systems using localized Lagrange multipliers[END_REF]. From a fundamental standpoint, the most satisfactory treatments have been based on variational principles, such as those cited in the Historical Background section. Such 'top down' derivations; however, may seem too abstract to engineers and scientists primarily concerned about how these methods can be useful in practical applications.

The aim of this paper is to deal with floating components in a more agile, step-by-step fashion by proceeding in a 'bottom up' manner. That is, starting with simple examples, and progressing toward the general case by stages while emphasizing the physical interpretation. Although linkage to variational principles is eventually made that is not an essential part of the exposition. It may be skipped by those readers primarily interested in modeling and implementation aspects.

Aside from the general goal of clarifying 'floating flexible dynamics', the paper has three specific objectives:

• Present the principal equations for a floating body or structure in a form suitable for matrix formulation and computer implementation. Those equations define that portion of the solution independent of constitutive properties and compatibility conditions. Mathematically that subspace spans the particular solution of the combined rigid motion and interface equilibrium equations. • Discuss the matrix form of the kinematic-static duality associated with the principal equations, which states that the self-equilibrium and rigid-body-basis matrices are the transposes of each other. • Apply the foregoing concepts to a set of examples of increasing generality.

Beyond the classical context of multibody, guidance and attitude dynamics, rigid-body motions and thus the d'Alembert-Lagrange principal equations play increasingly important roles in emerging fields such as molecular dynamics [START_REF] Frenkel | Understanding Molecular Simulation: From Algorithms to Applications[END_REF][START_REF] Konarev | MASSHA-a graphics system for rigid-body modelling of macromolecular complexes against solution scattering data[END_REF][START_REF] Nakano | A rigid-body-based multiple time scale molecular dynamics simulation of nanophase materials[END_REF][START_REF] Pirez | Molecular rigid-body displacements in a tetragonal lysozyme crystal confirmed by X-ray diffuse scattering[END_REF][START_REF] Sohlberg | Application of rigid-body dynamics and semiclassical mechanics to molecular bearings[END_REF], biophysics [START_REF] Dosset | A novel interactive tool for rigid-body modeling of multi-domain macromolecules using residual dipolar couplings[END_REF][START_REF] Knegtel | Molecular docking to ensembles of protein structures[END_REF][START_REF] May | Rigid body docking for virtual screening[END_REF][START_REF] Wang | Comparison of tRNA motions in the free and ribosomal bound structures[END_REF][START_REF] Wriggers | Protein domain movements: detection of rigid domains and visualization of hinges in comparisons of atomic coordinates[END_REF], as well as traditional applications such as multibody dynamics [START_REF] Fransen | Rigid body mode control of free-free launcher structural models in gravity fields[END_REF][START_REF] Loper | On viscous flow within a rotating spheroidal contained[END_REF][START_REF] O'brien | Synthesizing sounds from rigid-body simulations[END_REF][START_REF] Okuma | Identification of the principal rigid body modes under free-free boundary condition[END_REF].

The quasi-static case of the principal equations and equivalent forms have been utilized in the formulation of the classical force method (CFM) [START_REF] Argyris | Energy Theorems and Structural Analysis[END_REF][START_REF] Przemieniecki | Theory of Matrix Structural Analysis[END_REF], in the development of finite elements [START_REF] De Veubeke | Matrix structural analysis[END_REF][START_REF] Argyris | Die Methode der Finiten Elemente[END_REF], in the solution of singular systems for quasi-static problems [START_REF] Verchery | Régularisation du système de l' équilibre des structures élastiques discrètes[END_REF], in vibration eigenvalue problems [START_REF] Géradin | Mechanical Vibrations: Theory and Application to Structural Dynamics[END_REF], and for the parallel solution of large-scale structural systems [START_REF] Farhat | A method of finite-element tearing and interconnecting and its parallel solution algorithm[END_REF][START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF]. We believe that additional developments for multiphysics and multiscale applications remain to be done.

We emphasize that the focus of the present paper is on the d'Alembert-Lagrange principal equations (hereafter simply termed as principal equations), not on the more general d'Alembert-Lagrange equations or principles. The remainder of the paper is organized as follows.

Section 2 offers the historical background of the principal equations, from Newton to date, in the hope that the readers who are unfamiliar with the subject may want to start by reading this section while taking note of the references cited therein. It highlights roles played in the development of the force method of structural analysis, well-posedness of floating systems, subsystem interface conditions, and dual subdomain decomposition techniques for parallel computing.

In Section 3, we introduce the principal equations by way of a two-mass and three-spring planemotion problem. It is shown that the three principal equations are obtained by summing the forces and moments. We label the force and moment summation operators as self-equilibrium operators. Hence, the principal equations are obtained by the projection of the governing flexible equations of motion (EOM) onto the self-equilibrium operators. We show that the self-equilibrium and system rigid-body matrices are the transposes of each other. The section ends with a comment on an application of the self-equilibrium operators to establish solvability conditions for eigenproblems in which accelerations are prescribed.

Section 4 derives the self-equilibrium operator of a partially constrained static system, from which both the principal equations as well as conditions on admissible external forces readily emerge. Section 5 derives a general form of principal equations for partitioned linear flexible structural systems in conjunction with the method of localized Lagrange multipliers (LLM) [START_REF] Park | A variational framework for solution method developments in structural mechanics[END_REF][START_REF] Park | An algebraically partitioned FETI method for parallel structural analysis: algorithm description[END_REF][START_REF] Park | A variational principle for the formulation of partitioned structural systems[END_REF][START_REF] Park | A localized version of the method of Lagrange multipliers and its applications[END_REF]. Section 6 offers a step-by-step illustration of solving static partitioned floating systems by utilizing the principal equations.

Section 7 presents a discrete variational formulation of partitioned structural systems in terms of deformations, rigid-body modes (RBM), LLM, and interface frame displacements as independently varied state variables. The rigid-body mode equations are shown to conflate into the principal equations.

The Conclusions section lists potential future applications of the principal equations, with emphasis on the computational treatment of multiphysics and multiscale problems.

HISTORICAL BACKGROUND

This section outlines the original development of principal equations and some of their applications in classical and computational mechanics. No attempt to historical completeness is made. The reader may consult, for example, Dugas [START_REF] Dugas | A History of Mechanics[END_REF] or Truesdell and Toupin [START_REF] Truesdell | The classical field theories[END_REF] for classical developments or Farhat and Roux [START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF] for uses in domain decomposition methods.

Theoretical background

In their 1960 review of field theories Truesdell and Toupin [31, p. 595] state that 'The Principle of Virtual Work (PVW), which when the reaction of inertia is included may be called the Lagrange-d'Alembert principle, is the oldest general variational form of the equations of mechanics.' Lanczos [32, p. 92] remarks that 'd'Alembert's principle gives a complete solution of problems in mechanics. All different principles of mechanics are merely mathematically different formulations of d'Alembert's principle.' [Caution: the term 'principle' is not to be interpreted in the modern sense of a true variational principle such as Hamilton's. The PVW is actually what we call now a weak form or Galerkin form, with virtual displacements as test functions.]

As a variational statement d'Alembert's principle was first proven for a fully unconstrained system of bodies by Lagrange in his Mécanique Analytique [33, 2de Partie, 3me Section 3]. Thus a fair attribution should include both names. For example, in a Wikipedia article [34], it is called the d'Alembert-Lagrange principle, which is chronologically appropriate. D'Alembert's contribution to the dynamical PVW was indirect but important. He had asserted in his Trait é de Dynamique [START_REF]Traité de Dynamique[END_REF]Section 50] that an accelerating system can be transformed into an equivalent static system by adding the so-called inertial force and inertial torque or moment. The inertial force must act through the center of mass whereas the inertial torque can act anywhere. (The direct application of Newton's laws requires that the angular acceleration be applied only about the center of mass.) A coordinate system in which this substitution is valid is called an inertial frame. For example, if the total mass (assumed invariant) of the system is M and the center of mass acceleration with respect to an inertial frame is ẍ, the inertial force is -M ẍ.

Superposing the inertial force and torque reactions to the actual external forces and moments acting on the system gives the effective forces and effective moments acting on the body. (Some authors, such as Synge [START_REF] Synge | Classical dynamics[END_REF], reserve the name 'effective forces' for the inertial part.) Lagrange [33, 2de Partie, 1re Section] observes: 'Cette méthode réduit toutes les lois du mouvement des corps à celles de leur équilibre, et ramène ainsi la Dynamique à la Statique. ' An immediate consequence is the center-of-mass theorem: If a body is mass invariant and referred to an inertial frame then the center of mass moves as a mass-point, of mass equal to the mass of the body and subject to the effective force acting on the body [31, p. 532].T h i sl a w precedes the PVW by a century. It had been used for special cases by Newton [37, Corollary III to Laws of Motion], d'Alembert [START_REF]Traité de Dynamique[END_REF][START_REF]Traitéd el ' Équilibre et du Mouvement des Fluides pour servir de suite au Traitéd e Dynamique (2nd edn)[END_REF],a n dE u l e r [START_REF] Euler | Theoria Motus Corporum Solidorum seu Rigidrum ex Primis nostrae Cognitionis Principiis Stabilita et ad Omnis Motus, qui in hujusmodi Corpora Cadere Possunt[END_REF]. It was shown to apply to 'un corps quelconque' by Lagrange [33, 2de Partie, 3me Section 3].

Abetting the 'territorial split' of calculus dating back to Newton and Leibnitz, Euler, and Lagrange downplayed the use of Newton's laws. The equivalence of the PVW and Lagrange's analytical dynamics to Newton's laws and conservation principles was established nearly a century later by Thomson (Lord Kelvin) and Tait [START_REF] Thomson | Treatise on Natural Philosophy, Part I[END_REF]. Compared with Newton's laws, a major advantage of the PVW as an integral form of equations of balance is the automatic elimination of workless forces. In a continuum constrained body, the PVW does not include any stress required to maintain the constraints, as long as such stresses do not work. A classical example is the incompressibility constraint in isochoric motions. The PVW is, therefore, indeterminate with respect to workless forces or stress systems.

Substructuring, partitioning, and floating systems

Substructuring in the FEM context was introduced in the aerospace industry during the 1960s as a device to overcome resource limitations on computers of the time, as well as to stage communication between different design groups. Two flavors, based on the displacement and force methods, are well covered in [START_REF] Przemieniecki | Theory of Matrix Structural Analysis[END_REF]. For the force version, however, floating substructures were avoided apriori. This kind of spatial subdivision method evolved in three directions: superelements [START_REF] Egeland | SESAM-69: a general purpose finite element method program[END_REF] for multilevel decomposition in the late 1960s, partitions [START_REF] Park | Partitioned analysis of coupled systems[END_REF] for coupled problems in the late 1970s, and subdomains [START_REF] Farhat | A method of finite-element tearing and interconnecting and its parallel solution algorithm[END_REF][START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF] for parallel processing during the 1990s.

The treatment of floating bodies in the context of matrix structural analysis was systematically addressed by Fraeijs de Veubeke [START_REF] De Veubeke | Matrix structural analysis[END_REF]. Unfortunately this 1973 lecture material never appeared on a journal. Another pioneering presentation is that of Verchery [START_REF] Verchery | Régularisation du système de l' équilibre des structures élastiques discrètes[END_REF]. The treatment of interfaces by Lagrange multipliers as connectors has older roots. As a variational principle to treat actual discontinuities in continuum (field) models, it was first proposed by Prager [START_REF] Prager | Variational principles for linear elastostatics for discontinous displacements, strains and stresses[END_REF]. His principle links interface displacements through a single global Lagrange multiplier field. The idea was expanded and applied by Tong in the hybrid-FEM context [START_REF] Tong | New displacement finite element method for solid continua[END_REF]. It was systematically extended as a family of variational principles by Atluri [START_REF] Atluri | On 'hybrid' finite-element models in solid mechanics[END_REF], also with the development of hybrid FEM models in mind. Its use for parallel computation was successfully pursued by Farhat and Roux [START_REF] Farhat | A method of finite-element tearing and interconnecting and its parallel solution algorithm[END_REF][START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF], who used global Lagrange multipliers as subdomain connectors for structural and acoustic partitions in multiphysics applications. The use of LLM for interfacing substructures was proposed and developed in [START_REF] Park | A variational framework for solution method developments in structural mechanics[END_REF][START_REF] Park | An algebraically partitioned FETI method for parallel structural analysis: algorithm description[END_REF][START_REF] Park | A variational principle for the formulation of partitioned structural systems[END_REF][START_REF] Park | A localized version of the method of Lagrange multipliers and its applications[END_REF] and for fluid-structure interaction in [START_REF] Park | Partitioned formulation of internal fluid-structure interaction problems via localized Lagrange multipliers[END_REF][START_REF] Newton | Philosophiae Naturalis Principia Mathematica (Mathematical Principles of Natural Philosophy. Also known as the Principia)[END_REF].

Principal equations and the CFM

There is a link between the CFM of structural analysis and the principal equations that is worth noting. The separation of the internal force state into principal and deformational states has ancient roots in the CFM, tracing back to Maxwell, Betti, Castigliano, and Clapeyron in the XIX Century. The example of Section 3, which leads to Equations ( 4)-( 6), plainly illustrates that idea. The CFM decomposes the solution into a primary and a redundant-driven state. The former pertains to a fictitious, statically determinate 'base' structure brought into existence by physical cuts (partitions) that activate redundants. Adjusting the redundants to meet cut (partitioned) compatibility conditions yields the correct solution by superposition. It should, however, be pointed out that very few researchers, if any, have addressed the so-called principal state for the dynamical problems. As narrated in [START_REF] Felippa | A historical outline of matrix structural analysis: a play in three acts[END_REF], the chief technical difficulty that hampered the traditional approach in the computer era was the automated selection of redundancies. This obstacle has been attacked by several investigators since the mid-1960s. Notable efforts to date include the Gauss-Jordan rank technique [START_REF] Robinson | Structural Matrix Analysis for the Engineer[END_REF], minimum-norm solutions [START_REF] Kaneko | On computational procedures for the force method[END_REF][START_REF] Kaneko | Minimum norm solutions to linear elastic analysis problems[END_REF] suggested by the natural factor technique [START_REF] Argyris | The natural factor formulation of the stiffness matrix displacement method[END_REF],a n dt h e integrated force method [START_REF] Patnaik | An integrated force method for discrete analysis[END_REF][START_REF] Patnaik | The variational energy formulation for the integrated force method[END_REF].

A dynamic version of the force method

None of the foregoing force method variations has succeeded in the commercial FEM environment. There is a practical barrier: all of them require element reformulation in terms of flexibilities and element-by-element extraction of discrete equilibrium equations. That shift is unlikely to happen in view of the huge investment on stiffness-based element libraries as well as on element interoperability. There are also application range limitations. Unsupported mechanical systems do not possess conventional flexibilities. No dual method gracefully handles dynamic, stability, and kinematically non-linear problems, as primal methods effortlessly do. For example, treating inertial effects necessarily brings in flexibilities to pass from accelerations to forces as independent variables, entangling spurious elasticity properties into the realm of mass. We hope that the force decomposition method presented here in conjunction with LLM-as displayed, for example, in Equation (55)-forestalls those shortcomings and assists software reuse by retaining displacements and their time derivatives as primary variables.

A FLOATING SPRING-MASS SYSTEM

This section illustrates, through a simple example, the derivation of the principal equations for a floating system. Consider a spring-mass model of the boomerang shown in Figure 1, assumed to move in the x, y plane. [The model shown there is used only for illustrative purposes. For predicting actual motions it is totally inadequate as it lacks aerodynamic forces; furthermore, the beam stiffness of the boomerang is ignored, which leads to a spurious kinematic mode.]

The EOM are written first for a rigid system with instantaneous center of mass located at c, which is taken as origin of coordinates. Force balance in the x, y directions and momentum balance about the z direction give in which

m 0 0 m ü v + k x 0 0 k y u v = f x f y (1)
m = ⎡ ⎢ ⎣ m 1 00 0 m 2 0 00 m 3 ⎤ ⎥ ⎦ , k x = ⎡ ⎢ ⎣ k 1 -k 1 0 -k 1 k 1 0 00 0 ⎤ ⎥ ⎦ , k y = ⎡ ⎢ ⎣ 00 0 0 k 2 -k 2 0 -k 2 k 2 ⎤ ⎥ ⎦ u = ⎧ ⎪ ⎨ ⎪ ⎩ u 1 u 2 u 3 ⎫ ⎪ ⎬ ⎪ ⎭ , v = ⎧ ⎪ ⎨ ⎪ ⎩ v 1 v 2 v 3 ⎫ ⎪ ⎬ ⎪ ⎭ , f x = ⎧ ⎪ ⎨ ⎪ ⎩ f x 1 f x 2 f x 3 ⎫ ⎪ ⎬ ⎪ ⎭ , f y = ⎧ ⎪ ⎨ ⎪ ⎩ f y 1 f y 2 f y 3 ⎫ ⎪ ⎬ ⎪ ⎭ (2)

The principal equations

The foregoing EOM are recast in terms of d'Alembert's forces (f D ) by treating the inertial forces as if they are part of the applied forces:

f D = f D x f D y = f x -m ü-k x u f y -mv-k y v (3)
Translational and rotational balance can be further studied by collecting and adding up force and moment components. In the sequel matrix symbols S T x , S T y ,a n dS T will stand for the summation operators about x, y, and the rotation angle ≡ z about z, respectively. (The use of the transpose is for convenience in stating the static-kinematic duality discussed later.) For the (x, y) translational force balance we obtain

S T x f D x = f ext x -M üc = 0 with S T x =[111 ], M = (m 1 +m 2 +m 3 ) f ext x = ( f x 1 + f x 2 + f x 3 ), üc = (m 1 ü1 +m 2 ü2 +m 3 ü3 )/M (4)
S T y f D y = f ext y -M vc = 0 with S T y =[111 ], M = (m 1 +m 2 +m 3 ) f ext y = ( f y 1 + f y 2 + f y 3 ), vc = (m 1 v1 +m 2 v2 +m 3 v3 )/M (5)
For the angular moment balance about the mass center, denoting by ¨ , the rigid angular acceleration:

M D c = S T ⎧ ⎨ ⎩ f D x f D y ⎫ ⎬ ⎭ = M ext -J ¨ = 0 with S T =[-y 1c -y 2c -y 3c x 1c x 2c x 3c ], x ic = x i -x c , etc. M ext = S T f ext x f ext y , J = 3 i=1 m i (x 2 ic + y 2 ic ) üi v i = üc vc + -y ic x ic ¨ , 3 i=1 m i y ic = 0, 3 i=1 m i x ic = 0 (6)
Note that the summation operators cancel out the internal spring forces. This can be seen by introducing the internal force vector:

S T x f int x = 0, S T y f int y = 0, S T f int x f int y = 0 (7)
in which f int x = k x u and f int y = k y u are spring forces. Conditions ( 7) characterize the self-equilibrium states, also denoted as 'self-equilibrated states' in the literature. This terminology applies to floating structures or bodies without external or reaction forces. We will refer to (S T

x , S T y , S T ) as the selfequilibrium operators. The nontrivial overall balance equations for the example problem ( 4), ( 5), (6) will be referred to as d'Alembert-Lagrange principal equations, to recognize the pioneering contributions of d'Alembert [START_REF]Traité de Dynamique[END_REF] and Lagrange [START_REF] Lagrange | [END_REF] outlined in the Historical Background section. The abbreviation principal equations will be often used.

In FEM discretizations, freedoms are arranged node-by-node. This rearrangement can be done by the permutation matrix

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ u 1 v 1 u 2 v 2 u 3 v 3 ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ = P ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ u 1 u 2 u 3 v 1 v 2 v 3 ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 1 1 1 1 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ u 1 u 2 u 3 v 1 v 2 v 3 ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ (8)
where blank entries store zeros. Rearranging freedoms as per P and the d'Alembert forces as per f D , Equations ( 4)-( 6) can be compactly expressed as

S T f D = 0 in which S T = ⎡ ⎢ ⎢ ⎢ ⎣ S T x S T y S T ⎤ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎣ 101010 010101 -y 1c x 1c -y 2c x 2c -y 3c x 3c ⎤ ⎥ ⎥ ⎦ (9)
where d'Alembert's forces and the self-equilibrium operators are rearranged in accordance with the DOF permutations given in (8).

Physical interpretation

Plainly the rows of S T displayed in (9) are three independent, infinitesimal RBM of the example system. Thus, its transpose S is an RBM basis matrix as its columns span the space of rigid motions. Actual RBM amplitudes are obtained by using three dimensionless scaling coefficients:

x , y ,a n d rot , which are collected into the column 3-vector a. The actual rigid-body motion is given by u r = Sa, whereas the associated virtual rigid-body displacements are u r = S a. Applying the PVW to the example system yields

(f D ) T u r = u T r f D = a T S T f D = 0 ( 10
)
Since a is arbitrary, S T f D must vanish, which reproduces ( 9). The relationship is not accidental and can be shown to hold for any floating body, whether flexible or rigid:

The self-equilibrium and rigid-body-basis matrices are the transposes of each other. This is an instance of a static-kinematic duality: the equilibrium (balance) and kinematic operators are adjoints of each other (transposes in the case of real operators expressed as matrices). Note that the rigid motion basis collected in S spans, the infinitesimal RBMs and not the finite motions. This is evident from the PVW derivation: in the virtual displacement expression u r = S a, matrix S must be 'frozen' at the current configuration so variations in a merely scale that linearization.

The orthonormalized rigid-body matrix R is the column-orthogonalized S so that R T R = I.I n matrix form

R = SQ -1/2 with Q = S T S (11)
(In practice the Gram-Schmidt orthogonalization process or the Cholesky decomposition of Q may be used to produce R stably as long as S has full rank.) The RBM motion u r = Sa in terms of the R basis is u r = R a, with a = Q 1/2 a. For the example problem, the columns of S are already orthogonal if x 1c + x 2c + x 3c = 0a n dy 1c + y 2c + y 3c = 0, which holds if m 1 = m 2 = m 3 , and all that remains is normalizing columns of S to unit length. An alternative orthonormalization of R with respect to the mass matrix is often more convenient in dynamics.

Including deformations

The principal equations ( 9) model only rigid motions. As the example system has six DOF, three additional deformation modes must be included. The complete solution can be obtained by the superposition of rigid and deformation motions: u = u r +u d . These may be extracted from u using projection:

u d = P r u = (I-R R T )u, u r = u-u d = Sa = R a = R R T u (12)
Note that R T P r = P r R = 0. Furthermore P r = P T r and P 2 r = P r , showing that P r is an symmetric orthogonal projector. Additional mathematical properties of this decomposition are discussed in [START_REF] Felippa | The construction of free-free flexibility matrices as generalized stiffness inverses[END_REF][START_REF] Felippa | The construction of free-free flexibility matrices for multilevel structural analysis[END_REF]. Equation ( 12) may be interpreted by saying that the principal equations act as a kinematic constraint condition on the well-known flexible EOM, which for infinitesimal motions are

Ku+M ü = f (13)
For the example problem, the principal equations provide half of the total. For a large-scale problem, however, say, n≫10 5 , the principal equations provide at most six degrees of freedom, which is a tiny percentage. Thus, their constraint role can be appreciated. This role is further examined in Section 4 by considering two floating continuum substructures.

Quasi-static version and the vibration eigenproblem

For applications in which accelerations are given directly or indirectly it is convenient to recast [START_REF] May | Rigid body docking for virtual screening[END_REF] in the quasi-static format

Ku = f * with f * = f-M ü (14)
Here f * is called a statically equivalent force. If K is m-times singular with kernel (null space) spanned by the m columns of S (or R), Ku = f * is consistent if and only if

S T f * = 0 or R T f * = 0.
That is, the equivalent forces f * must be in self-equilibrium. If that solvability condition is not met, the right-hand side may be projected on the deformational motion space: f * d = P r f * , as explained by Fraeijs de Veubeke [START_REF] De Veubeke | Iteration in semidefinite eigenvalue problems[END_REF] and Geradin and Rixen [START_REF] Géradin | Mechanical Vibrations: Theory and Application to Structural Dynamics[END_REF].H o w e v e r ,Ku = f * d may still suffer from a singular K.I fR is known ap r i o r i , this can be overcome by using a stabilized stiffness:

K st u = P r (f d -M ü) in which K st = K+ R R T has full rank ( 15
)
Explicit formation of K st will result generally in a full matrix, hindering sparsity in large-scale problems. Sparsity may be retained, however, by using the techniques described in [START_REF] Felippa | The construction of free-free flexibility matrices for multilevel structural analysis[END_REF],i nw h i c h the original system ( 14) is constrained by 'penalty springs' adaptively injected during factorization, and solved for additional right-hand sides. The solution of ( 15) is formally u = Ff * ,i nw h i c hF = K -1 st P r = P r K -1 st is called the free-free flexibility. Technically, F is the Moore-Penrose generalized inverse of K. Because of the projector action u = Ff * will be simply u d .

An important application is the free-vibrations eigenproblem with f = 0 and üi = 2 i u i . Linear systems of the form Ku new = Mu old are repeatedly solved in Lanczos or subspace iteration eigensolvers. The stabilized form [START_REF] Wriggers | Protein domain movements: detection of rigid domains and visualization of hinges in comparisons of atomic coordinates[END_REF] is of interest should R, which collects eigenvectors for the m zero frequencies, be assembled ap r i o r iby geometric arguments, and spectral shifting not used.

A PARTLY UNCONSTRAINED STRUCTURE

To illustrate the principal equations for a partly unconstrained structure, consider the 4-bar mechanism linkage shown in Figure 2. (This is a generalization of an example from [56, p. 12],inwhich it is assumed that the bars are rigid.) Bars ( 1), (2), and (3) have identical lengths L, elastic moduli E, areas A, are pin-jointed at nodes 1-4, and constrained to move on the {x, y} plane. Inertial and damping effects are omitted for simplicity. The degrees of freedom are the {x, y} displacements of nodes 2 and 3, collected

in u =[u x2 u y2 u x3 u y3 ] T , with conjugate forces f =[ f x2 f y2 f x3 f y3 ] T .
Superscript 'ext' is used for given external forces. The arrays of bar axial forces and deformational elongations are denoted by p =[p (1) p (2) p (3) ] T and d =[d (1) d (2) d (3) ] T , respectively. We shall use the abbreviations c = cos and s = sin .If = 0 there is only one linearized rigid-body mode, which is the only column of R:

S = R =[-s -c -sc ] T (16)
This may be normalized as R = R/ √ 2. The 1×4 self-equilibrium matrix is S T = R T .

B = ⎡ ⎢ ⎣ c -s 00 -1 010 00 -c -s ⎤ ⎥ ⎦ , C = EA L I, K = EA L ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1+c 2 -sc -10
-sc s 2 00

-10 1 +c 2 sc 00sc s 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (17)
Observe that the static-kinematic duality also holds for deformational operators; this is not surprising as the problem is self-adjoint. Given that S T B T = 0 it follows that S T K = 0; whence the principal equation is

S T f D = S T (f ext -Ku) = S T f ext =-s( f ext x2 + f ext x3 )+c( f ext y3 -f ext y2 ) = 0
. This is a self-equilibrium condition that imposes a constraint on the external forces; e.g. any symmetric force system

f ext x2 =-f ext x3 , f ext y2 = f ext y3
satisfies it for any = 0. The inclusion of inertial and damping terms is straightforward. The case = 0, which has two independent linearized RBMs in 2D, requires special treatment and is omitted.

MULTIPLIER-CONNECTED FLOATING SUBSTRUCTURES

In this section we consider two substructures resulting from the partition of a semidiscrete, continuum finite element model undergoing 3D motions, as pictured in Figure 3 Here (i, j, k) are fixed unit base vectors attached to the inertial coordinate frame, whereas (u,v,w) are the displacement components. At the six nodal points

r i = d i +u i , d i = x i i+ y i j+ z i k, u i = u i i+v i j+w i k for i = 1, 2, 3,...,6 (19)
and similarly for ṙi , ri ,and r i . Kinematic node values are collected in vectors ū, ( u), ( ü) and ( ū), respectively, where the bar distinguishes the node displacements from the field in [START_REF] Okuma | Identification of the principal rigid body modes under free-free boundary condition[END_REF]. Volume forces are classified into applied body forces f E , internal forces f I due to elasticity, and inertia force r . Inter-partition interaction as well as ground reactions are treated with Lagrange multipliers . These are discretized into delta functions at boundary nodes and collected into a vector k.T h e configuration of this vector depends on details of the treatment of the interaction. Boundary freedoms affected by the interaction are collected in a vector u = B u in which B is a Boolean matrix that extracts u as B u. To account for possible interpartition dislocations and/or specified support motions the gap vector u g is introduced. Application of d'Alembert's principle for the partitioned system gives

= V ( r) T (f ext -f I -r) dV -{k T (B ū-u g )}=0 (20)
Substituting ( 18) and ( 19) into the above principle, carrying out the spatial integration, we obtain from the stationarity of the above expression:

ūT f D -k T {B ū-u g }=0, f D = fext u - fI u -B T k-m u ü fext u = V N T f E dV, fI u = V B T r dV, m u = V N T N dV e = B ū, u = 8 i=1 N i ( , , )u i = N ū ūT =[u 1 ,v 1 ,w 1 , u 2 ,v 2 ,w 2 ,...,u 8 ,v 8 ,w 8 ] (21)
where r and e are stresses and strains, respectively, derived from displacements, B is the standard strain-displacement relations in the finite element method, and N are displacement interpolation functions within the brick.

Principal equations

Adding up d'Alembert forces over all nodes gives

S T trans f D = S T trans { fext u -B T k-m u ü}=0
, where S T trans = I 1 I 2 ...

I 8 ( 22
)
where I is a (3×3) identity matrix and the fact that the internal forces cancel out for a floating structure, viz., S T trans fI u = 0, is utilized. The moment summation requires the following vectorial identities:

m z = r×f = ⎡ ⎢ ⎣ i j k ⎤ ⎥ ⎦ T ⎡ ⎢ ⎣ 0 -zy z 0 -x -yx 0 ⎤ ⎥ ⎦ ⎡ ⎢ ⎣ f x f y f z ⎤ ⎥ ⎦ , r = xi+ yj+ zk, f = f x i+ f y j+ f z k (23)
Let (x 0 , y 0 , z 0 ) denote the coordinates of the point about which moments are taken. Defining

[ di ] T = ⎡ ⎢ ⎣ 0 -(z i -z 0 )( y i -y 0 ) (z i -z 0 ) 0 -(x i -x 0 ) -(y i -y 0 )( x i -x 0 ) 0 ⎤ ⎥ ⎦ (24)
the moment summation can be compactly expressed as

S T rot f D = S T rot { fext u -B T k-m u ü}=0i n w h i c h S T rot =[ dT 1 dT 2 ... dT 8 ] (25)
The sum of moments of the internal forces also vanishes, viz., S T rot fI u = 0, regardless of the point about which of the moments are taken.

As discussed in the previous example, the principal Equations ( 22) and ( 25) model the admissible rigid-body motions and represent a constraint condition on the interface Lagrange multiplier k. This can be seen by writing

[S trans S rot ] T B T k =[S trans S rot ] T { fext u -m u ü} (26)
This gives

[S trans S rot ] T (B T k)-s = 0 whence s =[S trans S rot ] T { fext u -m u ü} (27)
Note that s = 0 represents the self-equilibrium equation if the system is completely free. Naturally, for a partitioned system that is a part of constrained system, it represents a violation of the constraint condition.

Moment summation for structures with rotational freedoms

For structures modeled with rotational degrees of freedom, a modification of the node-by-node force-moment summation is necessary. Assuming that nodal degrees of freedom at node i are arranged in the usual way:

(u i ,v i ,w i , xi , yi , zi ) then S T =[S trans S rot ] T = S T 1 S T 2 ... S T n S T i = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 100 0 0 0 010 0 0 0 001 0 0 0 0 -(z i -z 0 )( y i -y 0 ) 100
(z i -z 0 ) 0 -(x i -x 0 ) 010 -(y i -y 0 )( x i -x 0 ) 00 0 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (28)
Here (x 0 , y 0 , z 0 ) denotes a point around which all the moments are summed up for each given element or each partitioned substructure. These force-moment summation operators can be specialized to models that carry only selective translational and rotational degrees of freedom. For example, for a plane beam with three nodal degrees of freedom, (u,v, z ), the node-by-node summation operator is obtained by eliminating the rows and columns corresponding to (w, x , y ):

S T i = ⎡ ⎢ ⎣ 10 0 01 0 -(y i -y 0 )( x i -x 0 ) 1 ⎤ ⎥ ⎦ (29)
Other freedom configurations, such as used in plates and shell models, may be handled similarly.

Mass properties of arbitrary structures

A trivial yet practically important application of the d'Alembert-Lagrange principal equation is the computations of mass properties of complex geometry. To this end, we introduce the principal mass matrix, M , via the virtual work corresponding to the principal inertia force as

U I r = u T r m u ür , u r = Sa ⇓ U I r = a T M ä, M = S T m u S (30)
where we employed the summation operator (see Equations ( 22), [START_REF] Géradin | Mechanical Vibrations: Theory and Application to Structural Dynamics[END_REF], and ( 26)), and the mass matrix(m u ) is obtained by the standard FEM procedure. Specifically, M can be expressed as

M = M tt m tr m T tr J rr ( 31
)
where (M tt , J rr , m tr ) are the (3×3) translational mass, rotatory inertia, and coupling matrix between the translational and rotatory inertia matrix, which vanishes when the center of mass coincides with the point at which the moment self-equilibrium operator, S , is generated.

STAGED STATIC ANALYSIS

The following example introduces the idea of staged solution of a simple structural system subject to static loads and treated by Lagrange multiplier localization. The general formulation for the dynamic case is provided in the following section.

The LLM equations

Two fixed-free, parallel extensional springs of axial stiffness k 1 and k 2 that can move only along the x direction are pulled by given static forces as shown in Figure 4. Only x displacements are allowed and thus component subscripts are omitted. Each spring is taken to be a substructure, labeled (1) and ( 2). Substructures are localized with Lagrange multipliers as illustrated in Figure 4(b). As a result each spring 'floats'. Multipliers are connected to components of a force transducer device called a frame. Figure 4(c) shows the final FEM model with six displacements (four for the springs and two for the frame) and four multipliers. (The following section explains the notation used in Figure 4(c) in more detail.) The FEM model has four element-level stiffness equilibrium equations, four compatibility equations, and two multiplier equilibrium equations. These equations are organized, in that order, to produce the following matrix system, which is displayed both in compact block-matrix notation and in entry-level form:

⎡ ⎢ ⎢ ⎢ ⎣ KB 0 B T 0 -L 0 -L T 0 ⎤ ⎥ ⎥ ⎥ ⎦ ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ u k u f ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ k 1 -
⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ u (1) 1 u (1) 2 u (2) 1 u (2) 2 ... (1) 1
(1) 2

(2) 1

(2) 2 ...

u 2 ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ f (1) 1 =0 f (1) 2 =F 1 f (2) 1 =0 f (2) 2 =F 2 ........ 0 0 0 0 ........ f 2 = P ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ = ⎧ ⎪ ⎨ ⎪ ⎩ f 0 f f ⎫ ⎪ ⎬ ⎪ ⎭ (32)
In these matrix equations, the specified frame displacement u 1 = 0 has been explicitly removed. Substructure displacements are transformed to deformational-rigid normal coordinates using the modal relations denote deformational and rigid-body motion amplitudes, respectively, for substructure m. Those are collected in separate subvectors q and a. Applying the foregoing transformation to (32) produces the system

u = ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ u (1) 1 u (1) 2 u (2) 1 u (2) 2 ⎫ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎭ = ⎡ ⎢ ⎢ ⎢ ⎣ 10 | 10 -10| 10 01 | 01 0 -1 | 01 ⎤ ⎥ ⎥ ⎥ ⎦ ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ q (1) 1
q (2) 1 . . . (1) 1 (2) 1 ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ =[U R] q a (33)
⎡ ⎢ ⎢ ⎢ ⎢ ⎣ K d 0B d 0 00B r 0 B T d B T r 0 -L 00 -L T 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ q a k u f ⎫ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎭ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 2k
⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ q (1) 1 q (2) 1 ... (1) 1 (2) 1 ... (1) 1 (1) 2 (2) 1 (2) 2 ... u 2 ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ -F 1 -F 2 . . . F 1 F 2 . . . 0 0 0 0 . . . P ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ f d f r 0 f f ⎫ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎭ (34)
in which

K d = U T KU, B d = U T B, B r = U T R

Solution steps

The blocked system (34) will be solved in a hierarchical, step-by-step fashion. The static principal equations are obtained by grouping the second and fourth block-matrix equations, which involve only multipliers:

L k = B r -L T k = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 00 00 . . . . . . . . . . . . 0 -10-1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ (1) 1 (1) 2 (2) 1 (2) 2 ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ = ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ F 1 F 2 . . . P ¯⎫ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎭ = f r f f = f (35)
Denote by L + the pseudoinverse of L that satisfies the 4 Penrose conditions, and by Z ab a s i s for the null space of L ,t h a ti s ,L Z = 0. (If L is square and of full rank, L + = L -1 is the ordinary inverse while Z is void.) The general solution of ( 35) is the sum of the particular and homogeneous solutions: k = L + f +Z b,inwhichb is an array of parameters 1 ,..., k of length k equal to the null space dimension of L . For the example problem,

k = 1a n d k = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ (1) 1 (1) 2 (2) 1 (2) 2 ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ = L -G f +Z b = 1 2 ⎡ ⎢ ⎢ ⎢ ⎣ 3 1 - -1 3 1 - -1 ⎤ ⎥ ⎥ ⎥ ⎦ ⎧ ⎪ ⎨ ⎪ ⎩ F 1 F 2 P ⎫ ⎪ ⎬ ⎪ ⎭ + ⎡ ⎢ ⎢ ⎢ ⎣ 1 -1 -1 1 ⎤ ⎥ ⎥ ⎥ ⎦ 1 = 1 4 ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 3F 1 + F 2 +2P F 1 -F 2 -2P F 1 +3F 2 +2P -F 1 + F 2 -2P ⎫ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎭ + ⎡ ⎢ ⎢ ⎣ 1 -1 -1 1 ⎤ ⎥ ⎥ ⎦ 1 = k P + k (36)
We call k P the self-equilibrium principal solution.I fb is not void, the correction k = Z b is determined in the next step to account for flexibility and compatibility effects. By construction those two vectors are orthogonal: k T k P = 0, whence k P can be interpreted as the l 2 minimumnorm solution of the rigid motion and frame equilibrium equations. Thus, this first step has features in common with the minimum-norm approach to the force method [START_REF] Kaneko | On computational procedures for the force method[END_REF][START_REF] Kaneko | Minimum norm solutions to linear elastic analysis problems[END_REF]. The major difference is that no structural flexibilities are involved: as discussed in the following section B r and L are purely geometric matrices, and in fact L is Boolean.

Inserting k = k P + k in the first and third matrix equations and passing the terms involving k P to the RHS gives the symmetric system

K d B d Z Z T B T d 0 q b = ⎡ ⎢ ⎣ 2k 1 0 √ 2 02 k 2 - √ 2 √ 2 - √ 20 ⎤ ⎥ ⎦ ⎧ ⎪ ⎨ ⎪ ⎩ q (1) 1 q (2) 1 1 ⎫ ⎪ ⎬ ⎪ ⎭ =- 1 2 √ 2 ⎧ ⎪ ⎨ ⎪ ⎩ 3F 1 + F 2 +2P -F 1 +3F 2 +2P 0 ⎫ ⎪ ⎬ ⎪ ⎭ = f d -B d k P 0 (37)
Solving yields

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ q (1) 1 q (2) 1 1 ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ = 1 4(k 1 +k 2 ) ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ -2 √ 2(F 1 + F 2 + P) -2 √ 2(F 1 + F 2 + P) (k 1 -3k 2 )F 1 +(3k 1 -k 2 )F 2 +2(k 1 -k 2 )P ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ (38)
Replacing 1 in (36) furnishes the final form of the Lagrange multipliers

k = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ (1) 1 (1) 2 (2) 1 (2) 2 ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ = k P + k = 1 k 1 +k 2 ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ k 1 (F 1 + F 2 + P) k 2 F 1 -k 1 (F 2 + P) k 2 (F 1 + F 2 + P) -k 2 (F 1 + P)+k 1 F 2 ⎫ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎭ (39)
If computation of rigid body motion amplitudes and frame displacements is of no interest, as is occasionally the case in static analysis, the solution process may stop here. Else one proceeds to the third step, taking the compatibility equations: B d q+B r a-Lu f = 0, plugging in q from the second step, and solving for a and u f . Finally, q and a may be substituted in [START_REF] Lagrange | [END_REF] to recover physical displacements. For this example

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ (1) 1 (2) 1 u 2 ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ = 1 k 1 +k 2 ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ (F 1 + F 2 + P)/ √ 2 (F 1 + F 2 + P)/ √ 2
F 1 + F 2 + P ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ u (1) 1 u (1) 2 u (2) 1 u (2) 2 ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ = 1 √ 2 ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 10 1 0 -101 0 01 0 1 0 -101 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ q (1) 1 q (2) 1 (1) 1 (2) 1 ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ = ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ 0 u 2 0 u 2 ⎫ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎭ (40)
The foregoing three solution steps are flowcharted in Figure 5. It is highlighted there that the second step, which in large-scale problems accounts for the bulk of the computations, is readily parallelized since substructures, as well as supports treated as interfaces, are fully localized by the presence of the Lagrange multipliers.

The explicit reduction to substructural normal coordinates as displayed in Equation (34) was done here for exposition convenience. In actual large-scale implementations only the precalculation of the RBM basis R is necessary. Efficient preconditioned CG techniques may be used to solve (34) in parallel processing. The key observation is that once q is eliminated in (34) we can iteratively solve for k with the principal solution (35) as starting vector. In subsequent CG iterations, new iterates (as defined by b k at the kth iteration) are projected onto the principal solution basis L -G ; hence, subsequent iterates are orthogonal to the principal solution, i.e. they are purely deformational. Details of this solution procedure may be found in [START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF][START_REF] Park | An algebraically partitioned FETI method for parallel structural analysis: algorithm description[END_REF].

DISCRETE VARIATIONAL FORMULATION

This section extends the example-oriented derivations of the foregoing sections to a more general case: a linear elastic structure partitioned into N substructures, each discretized by FEM. Damping effects are ignored. Substructure meshes are assumed to match across interfaces. (The case of nonmatching meshes is more elaborate although it can be treated by similar methods; the reader is referred to [START_REF] Park | Interfacing nonmatching FEM meshes: the zero moment rule[END_REF][START_REF] Park | A simple algorithm for localized construction of nonmatching structural interfaces[END_REF][START_REF] Rebel | A contact formulation based on localised Lagrange multipliers: formulation and application to two-dimensional problems[END_REF] for details.)

To illustrate the notation introduced below, reference is often made to Figures 6 and7. The 2D finite element mesh of Figure 6(a) is divided into three substructures: 1 , 2 ,a n d 3 by the cuts shown in Figure 6(b). The supported (displacement prescribed) portion of the boundary is also cut off; hence, the three partitions are completely unconstrained. Nodes located along the cuts are called interface nodes. DOF allocated to interface nodes are called interface freedoms. Interface nodes are marked by ⊕ in Figure 6(c). Substructure nodes not located on interfaces are called internal; those are marked by in Figure 6(c).

A kinematic frame is inserted between substructure interfaces as well as between substructures and supports. This is illustrated in Figure 6(c), in which meshes, supports, and frame are drawn slightly offset as an aid to visualization. The frame is a purely geometric object placed over element sides (in 2D) or faces (in 3D). Unlike contact, fracture, or interface elements, it is not endowed with physical properties. For matching meshes, frame nodes are collocated with interfaced nodes.

The main purpose of the frame is to completely localize partition information, including interface forces. To this end, frame and interface nodes are linked with Lagrange multipliers as shown in Figure 6(c). Mathematically the multipliers are node-collocated delta functions whereas physically they are interface node forces. Choice of this discretization space circumvents the need to carry out integrals over interfaces.

Figure 7 shows global and local node identification. In the latter, nodes are labeled by an integer pair (i, n) that identifies the nth node of the ith partition; those are marked as i.n in Figure 7(b) for ease of reading. The frame is conventionally given index zero so (0, n) is simply shown as n there. When labeling node freedoms, the substructure index is placed as superscript for visualization convenience. Thus, u (3) x 7 denotes the x displacement of node 7 of substructure 3. Frame nodes may be labeled with subscript f if this helps to clarify the notation; for example, u f 7 denotes the subvector of node displacements of frame node 7. In superelement-like, multilevel versions of this approach (not considered here) a triplet (l, i, n) including a superelement level index l would be required for substructures, and a duplex for interface frames.

Localization

If all substructures are collectively assembled into a single structure, as in Figure 6(a), the overall node displacement vector, node velocity vector, node acceleration vector, node force vector, stiffness matrix, and mass matrix are denoted by u g , ug , üg , f g , K g ,andM g , respectively. The corresponding entities for substructure m are called u (m) , u(m) , ü(m) , f (m) , K (m) ,a n dM (m) for m = 1, 2,...,N , respectively. Symbols K and M, without superscripts, denote the block-diagonal matrices built up Multipliers are inserted for each interface node. Those pertaining to substructure m are collected in vector k (m) . This has the same length and arrangement as u exists. All multipliers are collected in vector k. The energy contribution of all multipliers can be compactly expressed as

d = k T (u b -B f u f ) (48)
This will be referred to as the interface energy potential. (Some authors, such as Fraeijs de Veubeke, call it a 'dislocation energy' following Friedrichs, who introduced the idea of 'constraint relaxation' in a different context.)

Equations of motion

The d'Alembert force vectors with damping ignored are f D g = f g -K g u g -M g üg for the assembly and

(f D i ) (m) = f (m) -K (m) u (m) -M (m) ü(m) (m = 1, 2,.
..,N ) for the mth substructure. The displacement-based discrete energy functional are denoted = (u g ) for the assembled structure and m = m (u m g ) for each individual substructure. The first variation of the global potential is g = u T g f D g . Setting this to zero yields the well-known global linear structural dynamics equation

M g üg +K g u g = f g (49)
For the partitioned system, the total energy is the sum of that contributed by individual substructures, plus the interface potential

(u g ) = N m=1 m (u m )+ N m=1 k m [(B m ) T u m -(B m f ) T u f ] = g (u g )+k m [B T u m -B T f u f ] (50)
In this functional u, u f ,andk are considered as independent variables and are varied independently. The first variation is

(u, u f , k) = u T f D + {k T (B T u-L f u f )} (51)
Setting the variation to zero yields the EOM ⎡

⎢ ⎢ ⎣ KB 0 B T 0 -L f 0 -L T f 0 ⎤ ⎥ ⎥ ⎦ ⎧ ⎪ ⎨ ⎪ ⎩ u k u f ⎫ ⎪ ⎬ ⎪ ⎭ + ⎡ ⎢ ⎣ M00 000 000 ⎤ ⎥ ⎦ ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ ü k ü f ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ = ⎧ ⎪ ⎨ ⎪ ⎩ f 0 0 ⎫ ⎪ ⎬ ⎪ ⎭ (52)
Note that K, which is a block-diagonal concatenation of singular stiffness matrices as per [START_REF] Egeland | SESAM-69: a general purpose finite element method program[END_REF], is highly rank deficient, since its kernel is spanned by the RBMs of all substructures. For example, if N = 1000 and each substructure has 6 RBMs, K will be 6000 times rank deficient. If the EOM (52) are solved by direct time integration, such deficiency is compensated by the presence of the mass matrix. But for certain situations described below, it is convenient to rewrite these equations explicitly in terms of deformational and rigid motions.

RBM-split EOM

The localized partitioned EOM (52) has full rank, provided the assembled total system Equation ( 49) has its full rank. However, this rank sufficiency is not adequate when one solves the equation in a partitioned manner. For example, when a subsystem becomes completely free-free as a result of partitioning, its self-equilibrium state (see, e.g. [START_REF] De Veubeke | Matrix structural analysis[END_REF]) must be part of the solution. The variational infusion of the self-equilibrium states for all the subsystems can be accomplished by decomposing the subsystem displacement into rigid and deformational: u = u r +u d . These are expressed in modal form u r = Ra and u d = Uq, where the modal basis matrices R and U are orthonormalized with respect to the mass matrix:

T M = M d , R T M = 0, KR = 0, R T MR = M r ( 53
)
in which if the eigenvectors are mass normalized, M d and M r become identity matrices. However, when the RBM are geometrically generated, they are not mass-normalized. Substitution of u = u r +u d = q+Ra into (51) gives

(q, a, u f , k ℓ ) = ( q+Ra) T f D ( q+Ra)+ [k T (B T ( q+Ra)-L f u f )] (54)
in which f D ( q+Ra) = f-K( q+Ra)-M( q+R ä). Setting = 0 and exploiting the orthonormality conditions (53) yields the rigid-deformational split EOM:

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ K d 0B d 0 00 B r 0 B T d B T r 0 -L f 00 -L T f 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ q a k u f ⎫ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎭ + ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ M d 00 0 0M r 00 00 0 0 00 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ q ä k ü f ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ = ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ f d f r 0 0 ⎫ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎭ (55) in which K d = U T KU, M d = U T MU, M r = R T MR, B d = U T B, B r = R T B, f d = U T f,a n df r = R T f.
The foregoing matrix equations have immediate physical interpretation: K d q+B d k+M d q = f d is the EOM for the substructural deformable modes, B r k+M r ä == f r governs substructural rigid motions, B T d q+B T r a-L f u f = 0 enforces the interface kinematic constraints, and finally L T f k = 0 expresses Newton's third law along the interfaces in terms of LLM. The first and second matrix equations can be localized to each substructure since they do not involve the frame displacements u f and consequently are fully parallelizable. The only linkage between the localized equations is provided by u f .

Remark

We have emphasized here the direct construction of self-equilibrium matrix operators (force and moment summation) as basis for building the principal equations. From these, RBM matrices follow by transposition. An alternative procedure is to try to construct the RBM matrices first, either by direct geometric arguments or by computer-based subdomain kinematic analysis based on the singular value decomposition (see, e.g. Farhat and Géradin [START_REF] Farhat | On the general solution by a direct method of a largescale singular system of linear equations: application to the analysis of floating structures[END_REF]). The latter technique has the advantage of accounting for the presence of internal mechanisms regardless of their origin (motion-control devices, spurious mesh modes, rank-deficient domain decomposition, etc.). If this alternative strategy is followed, the self-equilibrium conditions again follow from duality. The mass properties displayed in Equation ( 31) do not generally hold, however, for null space vectors that are extracted from singular stiffness matrices via an linear algebra procedure.

CONCLUSIONS

The focus of the paper has been the application of the d'Alembert-Lagrange principal equations to handle floating mechanical systems or components. As discussed in the Introduction, those modeling scenarios may arise directly in naturally unconstrained problems, such as a flying aircraft, or may occur as a result of a system decomposition process. In the latter case, the principal equations become a key ingredient of the formulation of the coupled discrete equations. We have shown that the process relies on PVW applied to the d'Alembert forces acting on the whole system, or on each individual component. Although the procedure can be followed as recipe by practicing engineers without a detailed knowledge of variational principles, work theorems or adjoint properties, such knowledge comes handy in understanding attributes such as the duality between static and kinematic operators.

Two important features of the principal equations presented here are their independence from constitutive and compatibility properties, and their indifference as to the source of forces acting on the system. The latter property represents an extension of the classical notion of 'self-equilibrium states' that refers to structures without external and reaction forces.

If the overall system is statically determinate, the principal equations directly supply the internal forces without the need of finding deformations and displacements. While in static analysis of fully supported structures this feature is common to flexibility-based methods, those cannot easily handle floating systems or dynamics. For statically indeterminate problems, the principal solution may serve as the initial step in an iterative process that injects compatibility and material properties in a hierarchical fashion. That process can be efficiently combined with Lagrange multiplier methods to address scenarios such as multiscale analysis or parallel computation.

As noted in the Historical Background section, multiplier-linked system decomposition methods have diverse roots in the past. In a taxonomy of model-based simulation for computational mechanics, the LLM methodology presented here may be categorized as a primal-dual hybrid: a combination of the displacement and force methods that tries to amalgam the best features of each while retaining application pliability. The following modeling and computational advantages can be cited for this combination.

1. If the system is undecomposed and fully supported, it reduces to the standard direct stiffness method, thus facilitating selective reuse of existing software. 2. At the other extreme, that of widespread cutting of rigid or highly stiff individual members it links naturally to multiplier-based multibody dynamics. 3. For structural models, it uses only displacement-connected elements in each substructure.

This protects the large existing investment of commercial FEM codes on element libraries, and retains interelement mixability as the tool to build complex models. 4. Cuts can be placed to evince interaction forces in problems where unexpected highly local effects call for deeper investigation (sort of an 'MRI structural test'). This is important in inverse problems for vulnerability/damage detection suggested by dynamic signal data. A related potential application is real-time, multisite hybrid testing, which is gaining importance as a result of networking advances.

5. Selective elimination of the frame reduces LLM to the mortar method for linking heterogeneous components, allowing result crosschecking with a developing body of multiphysics coupling technology. 6. As emphasized in the paper, the approach naturally handles floating systems and components.

Because displacements are an essential part of the formulation, it has no difficulties in modeling dynamics, vibration, and stability problems. 7. In problems that involve localized non-linearities (fracture, contact, etc.), cuts may be adaptively made to separate linear and non-linear portions of the overall system. 8. Lends itself to parallel implementation as well as treatment of nonmatching meshes. 9. It opens up a promising framework for partitioned analysis of coupled problems in multiphysics and multiscale modeling, as well as construction of reduced models by truncated representations of governing equations. Substantial developments in all of these areas, however, will be required to define correct and consistent interface conditions.

Two implementation obstacles remain as research areas that may be of interest to numerical analysts and applied mathematicians. The reliable automatic detection of non-RBM subsystem mechanisms (zero-energy internal modes) is notoriously difficult since it involves rank detection in floating-point arithmetic. Despite recent advances [START_REF] Farhat | On the general solution by a direct method of a largescale singular system of linear equations: application to the analysis of floating structures[END_REF], the topic deserves further study. Coupling subsystems or partitions of very different physical or inertial properties can lead to significant illconditioning in the interface equations. This can be expected to be particularly crucial in multiscale analysis. While ill-conditioning in structure-structure interaction has been alleviated with ad hoc scaling and filtering procedures [START_REF] González | Partitioned formulation of frictional contact problems using localized Lagrange multipliers[END_REF][START_REF] Park | On transforming tight coupling into loose coupling for a class of multiscale models[END_REF], a more general treatment seems appropriate.

Potential future applications of this approach include situations where noisy large motions, localization, contact-impact, heterogeneity, and/or scaling effects may disturb the smooth application of conventional displacement methods. Specific instances include multiscale inverse problems, space robotics, contact-impact, real-time testing, control of highly flexible systems, and obviously modeling and simulation of multiphysics and multiscale problems.

Figure 1 .

 1 Figure 1. A crude spring-mass 2D model of a boomerang. (Note this simple model does not capture boomerang motions.)

Figure 2 .

 2 Figure 2. A 4-bar plane mechanism linkage fabricated with three identical bars.

18 )Figure 3 .

 183 Figure 3. Two floating 3D substructures resulting from partition of a FEM model (the z-direction dimensions are suppressed for visual clarity): (a) assembled structure and (b) partition into two substructures.

Figure 4 .

 4 Figure 4. Two linear extensional springs in parallel, localized by Lagrange multipliers: (a) original structure; (b) localization; and (c) FEM model.

Figure 5 .

 5 Figure 5. Computational sequence followed in static analysis.

Figure 6 .

 6 Figure 6. Partition of a finite element mesh into 3 node-matching floating substructures, followed by introduction of interface frame and Lagrange multiplier connectors.

Figure 7 .

 7 Figure 7. Global and local node labeling.

  dot product (k (m) ) T u (m) b

  ................................................

	00 00 00k 2 -k 2 k 1 -k 1 k 1 00 -k 2 k 2	. . . 1000 . . . 0100 . . . 0010 . . . 0001	. . . 0 . . . 0 . . . 0 . . . 0
	1000	. . . 0000	. . . 0
	0100 0010	. . . 0000 . . . 0000	. . . -1 . . . 0
	0001	. . . 0000	. . . -1
	0000	. . . 0 -10-1	. . . 0

.................................................

  ............................................. ............................................. ................................................

	1 02 k 2 0 ...00	. . . 00 . . . 00 . . . 00	. . . . . . 00 -. . .	00 -00	. . . 0 . . . 0 . . . 0
	00 ...0	. . . 00 . . . 0	. . . 00 . . . 0000	. . . 0 . . . 0
	-0	0	. . . . . . 0	0	. . . 0000 . . . 0000	. . . -1 . . . 0
	0 -00	. . . 0 . . . 00	. . . 0000 . . . 0 -10-1	. . . -1 . . . 0

with the disconnected stiffness and mass matrices: (1) , K (2) ,...,K (N ) ], M = diag[M (1) , M (2) ,...,M (N ) ] (

Symbols u and f, without subscripts, denote the vector concatenation of the u (m) and f (m) , respectively. The substructure and global discretization levels are formally connected by Boolean matrices and are defined as

u (2) . . .

L (2) . . .

These can be compactly written as the supermatrix relations

The substructure-by-substructure assembly operation may be expressed as the congruences

although these entities are never actually assembled. For each substructure, we identify a subset u

Frame displacements are collected in vector u f . These will be considered as independent variables in the partitioned approach. Finally, we introduce N localization matrices B (m) f that link the frame displacements u f to the interface freedoms of the mth substructure as

Constraint treatment

The displacements of, say, frame node 2 and partition nodes 1.4 and 3.12 in Figure 7(b) must be the same. In a plane stress problem, this requires four constraints:

Note that multipliers are simply identified by the interface node label. The linkage to a specific frame node is fortuitous: it happens because we consider here only matching meshes.


Weiping Zhou

Arnaud Voituriez
email: arnaud.voituriez@cnrs.fr

Dr A Zhou

Voituriez

Asymmetric [3,3]-Sigmatropic Rearrangement of Sulfoniums for the Synthesis of Quaternary Center-Containing Polycyclic Molecules: Application to the Total Synthesis of Myrmenaphthol A

The asymmetric total synthesis of myrmenaphthol A, a natural product isolated from a Hawaiian sponge of the genus myrmekioderma, has been achieved in 12 steps. The key transformation of this synthesis is a gold-catalyzed [3,3]-sigmatropic rearrangement of sulfonium, starting from a chiral sulfoxide substrate and propargyl silane. The development of the methodology included the synthesis of benz[e]inden-2-one derivatives with quaternary centers (51-78% yield, 90-99% ee).

Steroids and their analogues are a unique class of natural products, found as active pharmaceutical ingredients (APIs) and present in many therapeutic agents. [1] Indeed, natural and synthetic steroids offer treatment possibilities for various cancers, inflammatory and cardiovascular diseases, among other pathologies. [2] In the present study, we were particularly interested in partially aromatic steroids (Scheme 1a). [3a] Undeniably, equilenin (6,8-didehydroestrone) (I) is the bestknown member of this subclass of natural molecules. In addition to being a component in conjugated estrogens, this product is the first sex hormone to be synthesized by organic chemists, in 1939 by the group of Bachmann. [3b] In medical chemistry, this tetracyclic skeleton has also proven to have very interesting biological properties, notably as estrogen receptor  agonist (ER EC50 = 21 nM). [4a] The molecule (II) has been synthesized by Micalizio and coworkers, via a metallacycle-mediated annulative cross-coupling reaction/cyclopropanation/rearrangement cascade process. [4] On the other hand, the natural product myrmenaphthol A 1, was isolated in 2019 from a Hawaiian sponge of the genus myrmekioderma. [5] This chiral molecule has some features that makes it a target of choice for organic chemists: (1) the presence of a quaternary stereocenter, (2) a hydroxyl group at the C2-position, while natural steroids, such as equinolin, estradiol or testosterone, have usually the C3 position oxidized;

(3) the presence of the naphthol core, which is unprecedented in a non-terrestrial natural steroidal product. Seduced by these features, we designed an asymmetric total synthesis of myrmenaphthol A 1 starting from readily available starting materials. In our retrosynthetic analysis (Scheme 1b), the free naphthol ring could be formed via demethylation of the corresponding methoxy group of (R)-7. The generation of the exocyclic olefin was designed by applying aldol condensation with acetaldehyde and subsequent dehydration. The enolic unit should be formed by an oxidation reaction of (R)-4f. As key step of this synthesis, we propose to synthesize (R)-4f using a new asymmetric gold(I)-catalyzed methodology, starting from proparyl silane and sulfoxide substrate (R)-2f. Finally, compound 2f would be synthesized in 3 steps from commercially available or easily accessible 6-methoxy-3,4-dihydrophenanthren-1(2H)-one S1. Scheme 1. Partially aromatic steroids, including myrmenaphthol A, and proposed strategies for its total synthesis.

The key reaction of this synthesis involves sulfonium intermediates and a [3,3]-sigmatropic rearrangement, [6] followed by an aldolization/elimination step (Scheme 1c). This approach is based on the precedents described in the literature on gold-catalyzed processes, [7] charge-accelerated [3,3]-rearrangements developed via interrupted Pummerer reaction/sigmatropic rearrangements [8] and our work in the synthesis of cyclopentenones with a C4-quaternary stereocenter. [9] The mechanism of this transformation includes several key-steps, such as: i) the gold-activation of propargyl silane; ii) the subsequent addition of vinyl sulfoxide 2 to the Au(I)activated alkyne leading to sulfonium intermediate (III); iii) [3,3]sigmatropic rearrangement to form the thionium ion and the hydrolysis and protodemetallation steps delivering the chiral compound 4' and finally iv) the aldol condensation leading to the desired products 4. As a preliminary work, we have envisaged first to apply this goldcatalyzed transformation to the formation of chiral tetrahydro-2Hcyclopenta[a]naphthalen-2-one derivatives 4a-e (Scheme 2). [START_REF] Nicolaou | Examples of quaternary center-containing [6.6.5]-tricarbocyclic skeleton in natural products: a)[END_REF][START_REF] Snyder | Asymmetric syntheses of such complex tricyclic molecules: a)[END_REF] Based on our previous results [9] and the Au(I)-catalyzed transformation depicted in Scheme 1c, the present study was initiated by using 10 mol% of (Johnphos)Au(MeCN)SbF6 gold catalyst and the enantioenriched sulfoxide substrate (S)-2a (99% ee), [START_REF] Mukaiyama | sulfur atom in the presence of mCPBA. b) Enantioenriched sulfoxide substrates: synthesis of vinyl bromide using Cu-mediated oxidation, starting from tetralones: b)[END_REF] in the presence of 2.0 equivalents of propargylic triisopropylsilane 3 at 30 °C. After 16 hours, the 1,4-diketone intermediate 4' was formed in more than 80% conversion by 1 H NMR. This seems to indicate a much higher reactivity of the cyclic olefinic bond in the [3,3]-sigmatropic rearrangement step, in comparison with precedent reactivity with acyclic vinylsulfoxides. [9b] Due to the complex mixture of 1,4-diketone containing both TIPS-substituted and desilylated products, the reaction mixture was dried under vacuum and mixed directly with NaOH in an ethanol/water solvent mixture, at 70 °C for 16 hours. The desired tricyclic compound (S)-4a was isolated in high yield (70% yield over 2 steps) and 90% ee (Scheme 2a). Increasing the reaction temperature to 40 °C gave similar results, in a shorter reaction time (16 h vs 36 h, 72% yield and 89% ee). Substrate (S)-2b, with an electron-donating methoxy group, exhibited a comparable reactivity, with a better chirality transfer (product 4b, 75% yield and 93% ee). Bulkier alkyl and aryl groups at the R 2 position, such as n-butyl, benzyl and phenyl substituents (substrates 3c-e) were well tolerated, providing the desired products 4c-e in 51-78% yields and 96-98% ees. It is interesting to note that substrates with longer alkyl chains or aryl groups gave a good reactivity, with a slightly better chirality transfer than vinyl sulfoxides with -methyl substituent (up to 98% ee vs 90% ee). Overall, the expected chiral polycyclic backbones 4a-e were synthesized in synthetically useful yields with high to excellent enantioselectivities. Importantly, the absolute configuration of vinylsulfoxides guides stereospecifically the direction of chirality transfer in the sulfonium intermediate and then determines the absolute configuration of quaternary centers in the ring-fused molecules 4.

[a] Reaction conditions: 2 (1.0 equiv), 3 (2.0 equiv), JohnphosAuSbF6•MeCN (10 mol %), water (3.0 equiv), dichloroethane (0.05 M) at 30-40 °C for 16-60 h. After solvent evaporation the residue was treated with NaOH (19 equiv) in a mixture EtOH/H2O, at 70 °C for 16 h. [b] Isolated yield after purification by column chromatography on silica gel. [c] Enantiomeric excesses determined by HPLC (High Performance Liquid Chromatography) analysis using a chiral stationary phase. Scheme 2. Substrate scope of the one-pot transformation for the formation of quaternary center-containing polycylic molecules 4a-f. [a] Consequently, starting with substrates (R)-2c and 2e gave after the one-pot transformation the corresponding tricyclic molecules (R)-4c and 4e in 52-72% yield and 95-99% ees (Scheme 2b). With the use of the optimized reaction conditions, substrate (R)-2f (98% ee) reacted at 40 °C to furnish after the rearrangement/aldol condensation sequence the key intermediate for the total synthesis of myrmenaphthol A. We had the great pleasure to isolate (R)-4f in 58% yield and 90% ee, at 0.10 mmol scale (Scheme 2c). It was again confirmed that the reaction temperature has a remarkable influence on the reaction rate (16 h at 40 °C vs 60 h at 30 °C). We then scaled up the synthesis of (R)-4f, using 0.32 mmol of vinyl sulfoxide (R)-2f and 7.5 mol% of gold complex. To our delight, a better isolated yield of 70% was obtained, while maintaining the same degree of chirality transfer. To continue the total synthesis of 1, and for the introduction of the enolic unit on the cyclopentenone part of the molecule, we first tried an epoxidation/opening reaction. Despite the different oxidants screened (H2O2, tBuO2H, DMDO) and in the presence of different base (NaOH, Triton B, K2CO3) at different temperatures, the desired product 7 was not obtained. We then decided to promote the epoxydation step by reducing the carbonyl function of 4f to the corresponding allylic alcohol 5 (Scheme 3). It is noteworthy that this alcohol is not stable on silica gel and it is directly transformed into the corresponding diene. Consequently, the crude allylic alcohol 5 was directly treated with m-CPBA, resulting in a 4:1 ratio of over-oxidized α,α'-dihydroxy ketone 6' and the epoxide ring opening α-hydroxy ketone 6. Due to the difficulty of separating these two compounds, the mixture was oxidized directly with Dess-Martin periodinane, to afford the formation of 1,2-dicarbonyl tertiary alcohol 7' and enol (R)-7 in the same 4:1 ratio. It was then envisaged that the deoxygenation of α-hydroxy ketone could promote the transformation of 7' to (R)-7. Fortunately, treatment with trimethylsilyl iodide (TMSI) [13] in chloroform converted the crude mixture to the expected -ketoenolic product (R)-7, in 54% yield over 4 steps. We subsequently envisaged to install the exocyclic double bond in compound 9, starting from 7 (Scheme 4). After deprotonation of enol (R)-7 with freshly prepared lithium diisopropylamide (LDA), dienolate intermediate was trapped by acetaldehyde, furnishing β-hydroxy ketone 8 in 61% yield as a 1/1.4 mixture of two diastereoisomers. The conversion was not complete and 35% of the starting keto-enolic compound 7 was recovered. Unfortunately, many attempts to improve this conversion have failed. The dehydration step was then performed under reflux with KOH in methanol, giving the Z/E mixture of the desired olefinic product 9 in a 1/3.8 ratio. Finally, it was decided to finish the total synthesis by direct demethylation of the methoxyl group, starting from this isomeric mixture. Unexpectedly, much effort has been devoted to the final demethylation step. No reaction occurred in the presence of BBr3, and partial or total decomposition was observed using BCl3 or TMSCl/NaI as reagents. However, a ray of hope has appeared with the use of pyridine•HCl at 165°C for 4 hours. Using such harsh reaction conditions allowed the isolation of myrmenaphthol A (E-1). The minor Z-isomer of the natural product was also isolated and characterized after separation on silica gel. 1 H and 13 C NMR spectra of the synthetic sample of myrmenaphthol A were consistent with those of the natural product. A negative optical rotation value of [α] D 20 -212 (c 0.40, CHCl3) was obtained, which is consistent with that of the literature (lit. [5a] : [α] D 20 -79 (c 0.2, CHCl3), despite a higher absolute value, confirming the (S) absolute configuration for the natural product.

To resume, Myrmenaphthol A 1 was obtained from 6-methoxy-3,4-dihydrophenanthren-1(2H)-one S1 in 12 chemical steps.

In conclusion, molecular complexity has been successfully introduced into the chemistry of asymmetric sulfonium [3,3]sigmatropic rearrangement, to reach quaternary stereocentercontaining complex polycyclic backbones. This methodology was validated by the synthesis of different tricyclic chiral compounds in 51-78% isolated yields and 90-99% ees. It is important to note the excellent chirality transfer from sulfur to carbon atom in all cases. The practicability of our new methodology was underlined by the first asymmetric total synthesis of myrmenaphthol A. In the future, it might be interesting and meaningful to examine this sulfonium chemistry with more complex molecules in order to establish an efficient toolbox to the synthesis of even more complex natural products.

A small library of tricyclic natural product-like molecules with quaternary stereocenters was rapidly synthesized. We observed an excellent chirality transfer from sulfur to carbon atom, furnishing the desired products in 90-99% enantiomeric excesses and 51-78% isolated yields. Furthermore, the feasibility of our new asymmetric methodology was highlighted by the asymmetric total synthesis of myrmenaphthol A

Scheme 3 .

 3 Scheme 3. Formation of the keto-enol compound (R)-7 from (R)-4f.

Acknowledgements

We thank the CSC (China Scholarship Council) for PhD grant to W. Z. and CNRS for financial support.

Keywords: Asymmetric • Sigmatropic rearrangement • benz [e]inden-2-one derivatives • Quaternary center • Myrmenaphthol A