Efficient Model Transition in Adaptive Multi-Resolution Modeling of Biopolymers

1. Introduction

Multibody dynamics methods are being used extensively to model biomolecular systems to study important physical phenomena occurring at different spatial and temporal scales [1], [2]. These systems may contain thousands or even millions of degrees of freedom, whereas, the size of the time step involved during the simulation is on the order of femto seconds. Examples of such problems may include proteins, DNAs, and RNAs. These highly complex physical systems are often studied at resolutions ranging from a fully atomistic model to coarse-grained molecules, up to a continuum level system [3], [4], [5]. In studying these problems, it is often desirable to change the system definition during the course of the simulation in order to achieve an optimal combination of accuracy and speed. For example, in order to study the overall conformational motion of a bimolecular system, a model based on super-atoms (beads) [6], [7] or articulated multi-rigid and/or flexible body [8], [9] can be used. Whereas, localized behavior has to be studied using fully atomistic models. In such cases, the need for the transition from a fine-scale to a coarse model and vice versa arises. Illustrations of a fully atomistic model of a molecule, and its coarse-grained model are shown in Fig. (1-a) and Fig. (1-b).

(a) Fully atomistic model	(b) Mixed type multibody model

Figure 1.

Illustration of a biomolecular system. a) Fully atomistic model. b) Coarse grain model with different rigid and flexible sub-domains connected to each other via kinematic joints.

Given the complexity and nonlinearity of challenging bimolecular systems, it is expected that different physical parameters such as dynamic boundary conditions and applied forces will have a significant affect on the behavior of the system. It is shown in [10] that time-invariant coarse models may provide inadequate or poor results and as such, an adaptive framework to model these systems should be considered [11]. Transitions between different system models can be achieved by intelligently removing or adding certain degrees of freedom $\left(dof\right)$. This change occurs instantaneously and as such, may be viewed as model changes as a result of impulsively applied constraints. For multi-rigid and flexible body systems, the transition from a higher fidelity (fine-scale) model to a lower fidelity model (coarse-scale) using divide-and-conquer algorithm (DCA) has been studied previously in [12], [13]. DCA efficiently provides the unique states of the system after this transition. In this chapter, we focus on the transition from a coarse model to a fine-scale model. When the system is modeled in an articulated multi-flexible-body framework, such transitions may be achieved by two different means. In the first, a fine-scale model is generated by adding flexible $dof$. This type of fine scaling may be necessary in order to capture higher frequency modes. For instance, when two molecules bind together, due to the impact, the higher frequency modes of the system are excited. The second type of fine scaling transition may be achieved through releasing the connecting joints in the multi-flexible-body system. In other words, certain constraints on joints are removed to introduce new $dof$ in the model.

In contrast to those types of dynamic systems in which the model definition is persistent, and the total energy of the system is conserved, the class of problems discussed here experiences discontinuous changes in the model definition and hence, the energy of the system must also change (nominally increase) in a discontinuous fashion. During the coarse graining process, based on a predefined metric, one may conclude that naturally existing higher modes are less relevant and can be ignored. As such, the kinetic energy associated with those modes must be estimated and properly accounted for, when transitioning back to the fine-scale model. Moreover, any change in the system model definition is assumed to occur as a result of impulsively applied constraints, without the influence of external loads. As such, the generalized momentum of the system must also be conserved [14]. In other words, the momentum of each differential element projected onto the space of admissible motions permitted by the more restrictive model (whether pre- or post-transition) when integrated over the entire system must be conserved across the model transition. If the generalized momentum is not conserved during the transition, the results are non-physical, and the new initial conditions for the rest of the simulation of the system are invalid.

In the next section, a brief overview of the DCA and analytical preliminaries necessary to the algorithm development are presented. The optimization problem associated with the coarse to fine scale transitioning is discussed next. Then the impulse-momentum formulation for transitioning from coarse models to fine-scale models in the articulated flexible-body scheme is presented. Finally conclusions are made.

2. Theoretical background

In this section, a brief introduction of the basic divide-and-conquer algorithm is presented. The DCA scheme has been developed for the simulation of general multi-rigid and multi-flexible-body systems [15], [16], [13], and systems with discontinuous changes [17], [12]. The basic algorithm described here is independent of the type of problem and is presented so that the chapter might be more self contained. In other words, it can be used to study the behavior of any rigid- and flexible-body system, even if the system undergoes a discontinuous change. Some mathematical preliminaries are also presented in this section which are important to the development of the algorithm.

2.1. Basic divide-and-conquer algorithm

The basic DCA scheme presented in this chapter works in a similar manner described in detail in [15], [16]. Consider two representative flexible bodies $k$ and $k+1$ connected to each other by a joint $J^k$ as shown in Fig. (2-a). The two points of interest, $H_1^k$ and $H_2^k$, on body $k$ are termed handles. A handle is any selected point through which a body interacts with the environment. In this chapter, we will limit our attention to each body having two handles, and each handle coincides with the joint location on the body, i.e. joint locations $J^{k-1}$ and $J^k$ in case of body $k$. Similarly, for body $k+1$, the points $H_1^{k+1}$ and $H_2^{k+1}$ are located at the joint locations $J^k$ and $J^{k+1}$, respectively. Furthermore, large rotations and translations in the flexible bodies are modeled as rigid body $dof$. Elastic deformations in the flexible bodies are modeled through the use of modal coordinates and admissible shape functions.

DCA is implemented using two main processes, hierarchic assembly and disassembly. The goal of the assembly process is to find the equations describing the dynamics of each body in the hierarchy at its two handles. This process begins at the level of individual bodies and adjacent bodies are assembled in a binary tree configuration. Using recursive formulations, this process couples the two-handle equations of successive bodies to find the two-handle equations of the resulting assembly. For example, body $k$ and body $k+1$ are coupled together to form the assembly shown in Fig. (2-b). At the end of the assembly process, the two-handle equations of the entire system are obtained.

(a) Consecutive bodies	(b) Assembly $k:k+1$

Figure 2.

Assembling of the two bodies to form a subassembly. a) Consecutive bodies k and k+1. b) A fictitious subassembly formed by coupling bodies k and k+1.

The hierarchic disassembly process begins with the solution of the two-handle equations associated with the primary node of the binary tree. The process works from this node to the individual sub-domain nodes of the binary tree to solve for the two-handle equations of the constituent subassemblies. This process is repeated until all unknowns (e.g., spatial constraint forces, spatial constraint impulses, spatial accelerations, jumps in the spatial velocities) of the bodies at the individual sub-domain level of the binary tree are known. The assembly and disassembly processes are illustrated in Fig. (?).

Figure 3.

The hierarchic assembly-disassembly process in DCA.

2.2. Analytical preliminaries

For convenience, the superscript $c$ shows that a quantity of interest is associated with the coarse model, while $f$ denotes that it is associated with the fine model. For example, the column matrix ${\left[\begin{array}{c}{v}_1\\ {\stackrel{�}{q}}_1\end{array}\right]}^c$ represents the velocity of handle-1 in the coarse model, and ${\left[\begin{array}{c}{v}_1\\ {\stackrel{�}{q}}_1\end{array}\right]}^f$ represents the velocity of the same handle in the fine-scale model. In these matrices, $v_1$ and ${\stackrel{�}{q}}_1$ are the spatial velocity vector of handle-1 and the associated generalized modal speeds, respectively.

As discussed previously, the change in the system model definition may occur by changing the number of flexible modes used to describe the behavior of flexible bodies, and/or the number of $dof$ of the connecting joints. To implement these changes in the system model mathematically, the joint free-motion map is defined as follows.

The joint free-motion map $P_R^{J^k}$ can be interpreted as the $6�{?}^k$ matrix of the free-modes of motion permitted by the ${?}^k$ degree-of-freedom joint, $J^k$. In other words, $P_R^{J^k}$ maps ${?}^k�1$ generalized speeds associated with relative free motion permitted by the joint into a $6�1$ spatial relative velocity vector which may occur across the joint, $J^k$ [15]. For instance, consider a transition in which a spherical joint in the system is altered, where only one $dof$ is locked about the first axis. The joint free-motion maps of the fine and coarse models in this case are shown in the following:

$$P_R^{J^{k}f}=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]\text{,}{P}_R^{J^{k}c}=\left[\begin{array}{cc}0& 0\\ 1& 0\\ 0& 1\\ 0& 0\\ 0& 0\\ 0& 0\end{array}\right]$$ uid10

We define the orthogonal complement of the joint free-motion map, $D_R^k$. As such, by definition one arises at the following

$${\left(D_R^k\right)}^T{P}_R^{J^k}=0$$ uid11

3. Optimization problem

Any violation in the conservation of the generalized momentum of the system in the transition between different models leads to non-physical results since the instantaneous switch in the system model definition is incurred without the influence of any external load. In other words, the momentum of each differential element projected onto the space of the admissible motions permitted by the more restrictive model (whether pre- or post-transition) when integrated over the entire system must be conserved across the model transition [14]. Jumps in the system partial velocities due to the sudden change in the model resolution result in the jumps in the generalized speeds corresponding to the new set of degrees of freedom. Since the model is instantaneously swapped, the position of the system does not change. Hence, the position dependent forces acting on the system do not change, and do not affect the generalized speeds. Any change in the applied loads (e.g., damping terms) which might occur due to the change in the model definition and the associated velocity jumps do not contribute to the impulse-momentum equations which describe the model transition. This is because these changes in the applied loads are bounded, and integrated over the infinitesimally short transition time.

Consider a fine-scale model with $n$ $dof$. Let the $dof$ of the model reduce to $n-m$ after the imposition of certain instantaneous constraints. In this case, the conservation of the generalized momentum of the system is expressed as

$$L^{c/c}=L^{f/c}$$ uid12

In the above equation, $L^{c/c}$ and $L^{f/c}$ represent the momenta of the coarse and fine models, respectively, projected on to the space of partial velocity vectors of the coarse model. Equation (?) provides a set of $n-m$ equations which are linear in the generalized speeds of the coarse model and solvable for the unique and physically meaningful states of the system after the transition to the coarser model.

Now consider the case in which, the coarse model is transitioned back to the fine-scale model. Equation (?) is still valid, and provides $n-m$ equations with $n$ unknown generalized speeds of the finer model. Furthermore, during the coarsening process, the level of the kinetic energy also drops because we chose to ignore certain modes of the system. However, in actual biomolecular systems such a decrease in energy does not happen. Consequently, it is important to realize the proper kinetic energy when transitioning back to the fine-scale model. Therefore, the following equation must be satisfied

$$K{E}^f=\frac{1}{2}{\left(u^f\right)}^T?u^f$$ uid13

In the above equation $u^f$ is the $n�1$ column matrix containing the generalized speeds of the fine model, and $?$ represents the generalized inertia matrix of the fine model. It is clear that Eqs.�(?) and (?) provide $n-m+1$ equations with $n$ unknowns. This indicates that the problem is under-determined when multiple $dof$ of the system are released. We may arrive at a unique or finite number of solutions, solving the following optimization problem

$$\begin{array}{ccc}& & \text{Optimize}J(u^f,t)\hfill \\ & & \text{Subjected}\text{to}{?}_i(u^f,t)=0,i=1,?,k\hfill \end{array}$$ uid14

In the above equation, $J$ is the physics- or knowledge- or mathematics-based objective function to be optimized (nominally minimized) subjected to the constraint equations ${?}_i$. In [18], [19], different objective functions are proposed for coarse to fine-scale transition problems. For instance, in order to prevent the generalized speeds of the new fine-scale model from deviating greatly from those of the coarse scale model, we may minimize the $L^2$ norm of the difference between the generalized speeds of the coarse and fine scale models as follows

$$J={(u^f-u^c)}^T(u^f-u^c)$$ uid15

As indicated previously, $(n-m)$ constraint equations governing the optimization problem are obtained from the conservation of the generalized momentum of the system within the transition. The rest of the constraint equations are obtained from other information about the system, such as the specific value of kinetic energy or the temperature of the system.

The generalized momenta balance equations from Eq.�(?) are expressed as

$$A{u}^f=b$$ uid16

where $A$ and $b$ are $(n-m)�n$ and $n$ known matrices, respectively, and $u^f$ is an $n�1$ column matrix of the generalized speeds of the fine-scale system model. As a part of the optimization problem, one must solve this linear system for $n-m$ dependent generalized speeds in terms of $m$ independent generalized speeds. Therefore, the optimization is performed on a much smaller number ($m$) variables, with a cost of $O\left(m^3\right)$. For a complex molecular system, $n$ could be very large, and $n>>m$, hence a significant reduction is achieved in the overall cost of optimization as compared to other traditional techniques, such as Lagrange multipliers [20]. However, the computations required to find the relations between dependent and independent generalized speeds can impose a significant burden on these simulations. It is shown in [21] that if traditional methods such as Gaussian elimination or LU factorization are used to find these relations, this cost tends to be $O\left(n{(n-m)}^2\right)$. The DCA scheme provided here finds these relations at the end of the hierarchic disassembly process with computational complexity of almost $O\left(n\right)$ in serial implementation. In other words, in this strategy, DCA formulates the impulse-momentum equations of the system which is followed by providing the relations between dependent and independent generalized speeds of the system in a timely manner. As such, this significantly reduces the costs associated with forming and solving the optimization problem in the transitions to the finer models.

4. DCA-based momenta balance for multi-flexible bodies

In this section, two-handle impulse-momentum equations of flexible bodies are derived. Mathematical modeling of the transition from a coarse model to a fine-scale model is discussed. For the fine-scale to coarse-scale model transition in multi-flexible-body system the reader is referred to [22], [23]. We will now derive the two-handle impulse-momentum equations when flexible degrees of freedom of a flexible body or the joints in the system are released. Then, the assembly of two consecutive bodies for which the connecting joint is unlocked is discussed. Finally, the hierarchic assembly-disassembly process for the multi-flexible-body system is presented.

4.1. Two-handle impulse-momentum equations in coarse to fine transitions

Now, we develop the two-handle impulse-momentum equations for consecutive flexible bodies in the transition from a coarse model to a fine-scale model. It is desired to develop the handle equations which express the spatial velocity vectors of the handles after the transition to the finer model as explicit functions of only newly introduced modal generalized speeds of the fine model. For this purpose, we start from the impulse-momentum equation of the flexible body as

$${?}^f{\mathrm{??}}_1^f-{?}^c{\mathrm{??}}_1^c={\left[\begin{array}{c}{?}_R\\ {?}_F\end{array}\right]}_1^c{?}_{t_c}^{t_f}{F}_{1c}dt+{\left[\begin{array}{c}{?}_R\\ {?}_F\end{array}\right]}_2^c{?}_{t_c}^{t_f}{F}_{2c}dt$$ uid18

where ${?}^f$ and ${?}^c$ are the inertia matrices associated with the fine-scale and coarse models, respectively. Also, $t_c$ and $t_f$ represent the time right before and right after the transition. The quantities ${?}_{t_c}^{t_f}{F}_{1c}dt$ and ${?}_{t_c}^{t_f}{F}_{2c}dt$ are the spatial impulsive constraint forces on handle-1 and handle-2 of the flexible body. The matrices ${\left[\begin{array}{c}{?}_R\\ {?}_F\end{array}\right]}_1^c$ and ${\left[\begin{array}{c}{?}_R\\ {?}_F\end{array}\right]}_2^c$ are the coefficients resulting from the generalized constraint force contribution at handle-1 and handle-2, respectively. Moreover, in Eq.�(?), the impulses due to the applied loads are not considered since they represent a bounded loads integrated over an infinitesimal time interval. For detailed derivation of these quantities the reader is referred to [13]. It is desired to develop the handle equations which provide the spatial velocity vectors of the handles right after the transition to the fine-scale model in terms of newly added modal generalized speeds. Therefore, in Eq.�(?), the inertia matrix of the flexible body is represented by its components corresponding to rigid and flexible modes, as well as the coupling terms

$$\begin{array}{ccc}\hfill {\left[\begin{array}{cc}{?}_{RR}& {?}_{RF}\\ {?}_{FR}& {?}_{FF}\end{array}\right]}^f{\left[\begin{array}{c}{v}_1\\ \stackrel{�}{q}\end{array}\right]}^f& =& {\left[\begin{array}{c}{?}_R\\ {?}_F\end{array}\right]}_1^c{?}_{t_c}^{t_f}{F}_{1c}dt+{\left[\begin{array}{c}{?}_R\\ {?}_F\end{array}\right]}_2^c{?}_{t_c}^{t_f}{F}_{2c}dt\hfill \\ & +& {\left[\begin{array}{cc}{?}_{RR}& {?}_{RF}\\ {?}_{FR}& {?}_{FF}\end{array}\right]}^c{\left[\begin{array}{c}{v}_1\\ \stackrel{�}{q}\end{array}\right]}^c\hfill \end{array}$$ uid19

which is decomposed to the following relations

$$\begin{array}{ccc}\hfill {?}_{FF}^f{\stackrel{�}{q}}^f& =& {?}_{F1}^c{?}_{t_c}^{t_f}{F}_{1c}dt+{?}_{F2}^c{?}_{t_c}^{t_f}{F}_{2c}dt+{?}_{FR}^c{v}_1^c+{?}_{FF}^c{\stackrel{�}{q}}^c-{?}_{FR}^f{v}_1^f\hfill \\ \hfill {?}_{RR}^f{v}_1^f& =& {?}_{R1}^c{?}_{t_c}^{t_f}{F}_{1c}dt+{?}_{R2}^c{?}_{t_c}^{t_f}{F}_{2c}dt+{?}_{RR}^c{v}_1^c+{?}_{RF}^c{\stackrel{�}{q}}^c-{?}_{RF}^f{\stackrel{�}{q}}^f\hfill \end{array}$$ uid20

Since the generalized momentum equations are calculated based on the projection onto the space of the coarser model, the matrix ${?}^f$ is not a square matrix and thus ${?}_{FF}^f$ is not invertible. However, we can partition Eq.�(?) in terms of dependent (those associated with the coarser model) and independent (newly introduced) generalized speeds as

$$\begin{array}{ccc}\hfill \left[\begin{array}{ccc}{?}_{FF}^{f_d}& ?& {?}_{FF}^{f_i}\end{array}\right]\left[\begin{array}{c}{\stackrel{�}{q}}^{f_d}\\ ?\\ {\stackrel{�}{q}}^{f_i}\end{array}\right]& =& {?}_{F1}^c{?}_{t_c}^{t_f}{F}_{1c}dt+{?}_{F2}^c{?}_{t_c}^{t_f}{F}_{2c}dt\hfill \\ & +& {?}_{FR}^c{v}_1^c+{?}_{FF}^c{\stackrel{�}{q}}^c-{?}_{FR}^f{v}_1^f\hfill \end{array}$$ uid21

Using the above relation, the expression for the dependent generalized modal speeds is written as

$$\begin{array}{ccc}\hfill {\stackrel{�}{q}}^{f_d}& =& {\left({?}_{FF}^{f_d}\right)}^{-1}[{?}_{F1}^c{?}_{t_c}^{t_f}{F}_{1c}dt+{?}_{F2}^c{?}_{t_c}^{t_f}{F}_{2c}dt+{?}_{FR}^c{v}_1^c\hfill \\ & +& {?}_{FF}^c{\stackrel{�}{q}}^c-{?}_{FR}^f{v}_1^f-{?}_{FF}^{f_i}{\stackrel{�}{q}}^{f_i}]\hfill \end{array}$$ uid22

Defining

$$\begin{array}{ccc}\hfill {?}_{RF}^f& =& \left[\begin{array}{ccc}{?}_{RF}^{f_d}& ?& {?}_{RF}^{f_i}\end{array}\right]\hfill \\ \hfill ?& =& {[{?}_{RR}^f-{?}_{RF}^{f_d}{\left({?}_{FF}^{f_d}\right)}^{-1}{?}_{FR}^f]}^{-1}\hfill \\ \hfill {?}_1& =& [{?}_{R1}^c-{?}_{RF}^{f_d}{\left({?}_{FF}^{f_d}\right)}^{-1}{?}_{F1}^c]\hfill \\ \hfill {?}_2& =& [{?}_{R2}^c-{?}_{RF}^{f_d}{\left({?}_{FF}^{f_d}\right)}^{-1}{?}_{F2}^c]\hfill \\ \hfill {?}_3& =& [{?}_{RR}^c-{?}_{RF}^{f_d}{\left({?}_{FF}^{f_d}\right)}^{-1}{?}_{FR}^c]v_1^c+[{?}_{RF}^c-{?}_{RF}^{f_d}{\left({?}_{FF}^{f_d}\right)}^{-1}{?}_{FF}^c]{\stackrel{�}{q}}^c\hfill \\ \hfill {?}_4& =& [{?}_{RF}^{f_d}{\left({?}_{FF}^{f_d}\right)}^{-1}{?}_{FF}^{f_i}-{?}_{RF}^{f_i}]\hfill \\ \hfill {?}_{1i}& =& ?{?}_i,(i=1,2,3,4)\hfill \end{array}$$ uid23

and using Eqs.�(?) and (?), the spatial velocity vector of handle-1 can be written in terms of the independent modal speeds

$$v_1^f={?}_{11}{?}_{t_c}^{t_f}{F}_{1c}dt+{?}_{12}{?}_{t_c}^{t_f}{F}_{2c}dt+{?}_{13}+{?}_{14}{\stackrel{�}{q}}^{f_i}$$ uid24

As such, the spatial velocity vector of handle-2 becomes

$$v_2^f={\left(S^{k1k2}\right)}^T{v}_1^f+{?}_2^f{\stackrel{�}{q}}^f$$ uid25

Employing the same partitioning technique, Eqs.�(?) can be written as

$$\begin{array}{ccc}\hfill v_2^f& =& {\left(S^{k1k2}\right)}^T{v}_1^f+\left[\begin{array}{ccc}{?}_2^{f_d}& ?& {?}_2^{f_i}\end{array}\right]\left[\begin{array}{c}{\stackrel{�}{q}}^{f_d}\\ ?\\ {\stackrel{�}{q}}^{f_i}\end{array}\right]\hfill \\ \hfill ?v_2^f& =& {\left(S^{k1k2}\right)}^T{v}_1^f+{?}_2^{f_d}{\stackrel{�}{q}}^{f_d}+{?}_2^{f_i}{\stackrel{�}{q}}^{f_i}\hfill \end{array}$$ uid26

Using

$$\begin{array}{ccc}\hfill {?}_{21}& =& [{\left(S^{k1k2}\right)}^T{?}_{11}+{?}_2^{f_d}{\left({?}_{FF}^{f_d}\right)}^{-1}{?}_{F1}^c-{?}_2^{f_d}{\left({?}_{FF}^{f_d}\right)}^{-1}{?}_{FR}^f{?}_{11}]\hfill \\ \hfill {?}_{22}& =& [{\left(S^{k1k2}\right)}^T{?}_{12}+{?}_2^{f_d}{\left({?}_{FF}^{f_d}\right)}^{-1}{?}_{F2}^c-{?}_2^{f_d}{\left({?}_{FF}^{f_d}\right)}^{-1}{?}_{FR}^f{?}_{12}]\hfill \\ \hfill {?}_{23}& =& {\left[S^{k1k2}\right)}^T{?}_{13}+{?}_2^{f_d}{\left({?}_{FF}^{f_d}\right)}^{-1}{?}_{FR}^c{v}_1^c+{?}_2^{f_d}{\left({?}_{FF}^{f_d}\right)}^{-1}{?}_{FF}^c{\stackrel{�}{q}}^c\hfill \\ & -& {?}_2^{f_d}{\left({?}_{FF}^{f_d}\right)}^{-1}{?}_{FR}^f{?}_{13}]\hfill \\ \hfill {?}_{24}& =& {\left[S^{k1k2}\right)}^T{?}_{14}+{?}_2^{f_i}-{?}_2^{f_d}{\left({?}_{FF}^{f_d}\right)}^{-1}{?}_{FF}^{f_i}-{?}_2^{f_d}{\left({?}_{FF}^{f_d}\right)}^{-1}{?}_{FR}^f{?}_{14}]\hfill \end{array}$$ uid27

and Eq.�(?), the spatial velocity vector of handle-2 can be written as

$$v_2^f={?}_{21}{?}_{t_c}^{t_f}{F}_{1c}dt+{?}_{22}{?}_{t_c}^{t_f}{F}_{2c}dt+{?}_{23}+{?}_{24}{\stackrel{�}{q}}^{f_i}$$ uid28

Equations (?) and (?) are now in two-handle impulse-momentum form and along with Eq.�(?), give the new velocities associated with each handle after the transition. These equations express the spatial velocity vectors of the handles of the body as well as the modal generalized speeds which have not changed within the transition in terms of the newly added modal generalized speeds. This important property will be used in the optimization problem to provide the states of the system after the transition to the finer models.

As such, the two-handle equations describing the impulse-momentum of two consecutive bodies, body $k$ and body $k+1$ are expressed as

$$\begin{array}{ccc}\hfill v_1^{\left(k\right)f}& =& {?}_{11}^{\left(k\right)}{?}_{t_c}^{t_f}{F}_{1c}^{\left(k\right)}dt+{?}_{12}^{\left(k\right)}{?}_{t_c}^{t_f}{F}_{2c}^{\left(k\right)}dt+{?}_{13}^{\left(k\right)}+{?}_{14}^{\left(k\right)}{\stackrel{�}{q}}^{\left(k\right)f_i}\hfill \\ \hfill v_2^{\left(k\right)f}& =& {?}_{21}^{\left(k\right)}{?}_{t_c}^{t_f}{F}_{1c}^{\left(k\right)}dt+{?}_{22}^{\left(k\right)}{?}_{t_c}^{t_f}{F}_{2c}^{\left(k\right)}dt+{?}_{23}^{\left(k\right)}+{?}_{24}^{\left(k\right)}{\stackrel{�}{q}}^{\left(k\right)f_i}\hfill \\ \hfill v_1^{(k+1)f}& =& {?}_{11}^{(k+1)}{?}_{t_c}^{t_f}{F}_{1c}^{(k+1)}dt+{?}_{12}^{(k+1)}{?}_{t_c}^{t_f}{F}_{2c}^{(k+1)}dt+{?}_{13}^{(k+1)}+{?}_{14}^{(k+1)}{\stackrel{�}{q}}^{(k+1)f_i}\hfill \\ \hfill v_2^{(k+1)f}& =& {?}_{21}^{(k+1)}{?}_{t_c}^{t_f}{F}_{1c}^{(k+1)}dt+{?}_{22}^{(k+1)}{?}_{t_c}^{t_f}{F}_{2c}^{(k+1)}dt+{?}_{23}^{(k+1)}+{?}_{24}^{(k+1)}{\stackrel{�}{q}}^{(k+1)f_i}\hfill \end{array}$$ uid29

4.2. Assembly process and releasing the joint between two consecutive bodies

In this section, a method to combine the two-handle equations of individual flexible bodies to form the equations of the resulting assembly is presented. Herein, the assembly process of the consecutive bodies is discussed only within the transition from a coarse model to a finer model. This transition is achieved by releasing the joint between two consecutive bodies. Clearly, this would mean a change in the joint free-motion map $P_R^{J^k}$ and its orthogonal complement $D_R^{J^k}$. It will become evident that the assembly process of the consecutive bodies for the fine to coarse transition is similar, and the associated equations can be easily derived by following the given procedure.

From the definition of joint free-motion map, the relative spatial velocity vector at the joint between two consecutive bodies is expressed by the following kinematic constraint

$$v_1^{(k+1)f}-v_2^{\left(k\right)f}=P_R^{J^{k}f}{u}^{(k/k+1)f}$$ uid31

In the above equation, $u^{(k/k+1)f}$ is the relative generalized speed defined at the joint of the fine model. From Newton's third law of motion, the impulses at the intermediate joint are related by

$${?}_{t_c}^{t_f}{F}_{2c}^{\left(k\right)}dt=-{?}_{t_c}^{t_f}{F}_{1c}^{(k+1)}dt$$ uid32

Substituting Eqs.�(?), (?) and (?) into Eq.�(?) results in

$$\begin{array}{ccc}\hfill ({?}_{11}^{(k+1)}& +& {?}_{22}^{\left(k\right)}){?}_{t_c}^{t_f}{F}_{1c}^{(k+1)}dt={?}_{21}^{\left(k\right)}{?}_{t_c}^{t_f}{F}_{1c}^{\left(k\right)}dt-{?}_{12}^{(k+1)}{?}_{t_c}^{t_f}{F}_{2c}^{(k+1)}dt\hfill \\ & +& {?}_{23}^{\left(k\right)}-{?}_{13}^{(k+1)}+{?}_{24}^{\left(k\right)}{\stackrel{�}{q}}^{\left(k\right)f_i}-{?}_{14}^{(k+1)}{\stackrel{�}{q}}^{(k+1)f_i}+P_R^{J^{k}f}{u}^{(k/k+1)f}\hfill \end{array}$$ uid33

Using the definition of the joint free-motion map, the spatial constraint impulses lie exactly in the space spanned by the orthogonal complement of joint free-motion map of the coarser model. These constraint impulses can be expressed as

$${?}_{t_c}^{t_f}{F}_{1c}^{(k+1)}dt=D_R^{J^{k}c}{?}_{t_c}^{t_f}{\mathrm{??}}_{1c}^{(k+1)}dt$$ uid34

In the above equation, ${?}_{t_c}^{t_f}{\mathrm{??}}_{1c}^{(k+1)}dt$ is an ordered measure number of the impulsive constraint torques and forces. Pre-multiplying Eq.�(?) by ${\left(D_R^{J^{k}c}\right)}^T$, one arrives at the expression for ${?}_{t_c}^{t_f}{F}_{1c}^{(k+1)}dt$ as

$$\begin{array}{ccc}\hfill {?}_{t_c}^{t_f}{F}_{1c}^{(k+1)}dt& =& X{?}_{21}^{\left(k\right)}{?}_{t_c}^{t_f}{F}_{1c}^{\left(k\right)}dt-X{?}_{12}^{(k+1)}{?}_{t_c}^{t_f}{F}_{2c}^{(k+1)}dt\hfill \\ & +& XY+X{?}_{24}^{\left(k\right)}{\stackrel{�}{q}}^{\left(k\right)f_i}-X{?}_{14}^{(k+1)}{\stackrel{�}{q}}^{(k+1)f_i}+X{P}_R^{J^{k}f}{u}^{(k/k+1)f}\hfill \end{array}$$ uid35

where

$$\begin{array}{ccc}\hfill X& =& D_R^{J^{k}c}{\left[{\left(D_R^{J^{k}c}\right)}^T({?}_{11}^{(k+1)}+{?}_{22}^{\left(k\right)})D_R^{J^{k}c}\right]}^{-1}{\left(D_R^{J^{k}c}\right)}^T\hfill \\ \hfill Y& =& {?}_{23}^{\left(k\right)}-{?}_{13}^{(k+1)}\hfill \end{array}$$ uid36

Using Eqs.�(?), (?), and (?), we write the two-handle equations for the assembly $k:k+1$

$$\begin{array}{ccc}\hfill v_1^{(k:k+1)f}& =& {?}_{11}^{(k:k+1)}{?}_{t_c}^{t_f}{F}_{1c}^{\left(k\right)}dt+{?}_{12}^{(k:k+1)}{?}_{t_c}^{t_f}{F}_{2c}^{(k+1)}dt\hfill \\ & +& {?}_{13}^{(k:k+1)}+{?}_{14}^{(k:k+1)}{\stackrel{�}{q}}^{\left(k\right)f_i}+{?}_{15}^{(k:k+1)}{\stackrel{�}{q}}^{(k+1)f_i}+{?}_{16}^{(k:k+1)}{u}^{(k/k+1)f}\hfill \\ \hfill v_2^{(k:k+1)f}& =& {?}_{21}^{(k:k+1)}{?}_{t_c}^{t_f}{F}_{1c}^{\left(k\right)}dt+{?}_{22}^{(k:k+1)}{?}_{t_c}^{t_f}{F}_{2c}^{(k+1)}dt\hfill \\ & +& {?}_{23}^{(k:k+1)}+{?}_{24}^{(k:k+1)}{\stackrel{�}{q}}^{\left(k\right)f_i}+{?}_{25}^{(k:k+1)}{\stackrel{�}{q}}^{(k+1)f_i}+{?}_{26}^{(k:k+1)}{u}^{(k/k+1)f}\hfill \end{array}$$ uid37

where:

$$\begin{array}{ccc}\hfill {?}_{11}^{(k:k+1)}& =& {?}_{11}^{\left(k\right)}-{?}_{12}^{\left(k\right)}X{?}_{21}^{\left(k\right)}\hfill \\ \hfill {?}_{12}^{(k:k+1)}& =& {?}_{12}^{\left(k\right)}X{?}_{12}^{(k+1)}\hfill \\ \hfill {?}_{13}^{(k:k+1)}& =& {?}_{13}^{\left(k\right)}-{?}_{12}^{\left(k\right)}XY\hfill \\ \hfill {?}_{14}^{(k:k+1)}& =& {?}_{14}^{\left(k\right)}-{?}_{12}^{\left(k\right)}X{?}_{24}^{\left(k\right)}\hfill \\ \hfill {?}_{15}^{(k:k+1)}& =& {?}_{12}^{\left(k\right)}X{?}_{14}^{(k+1)}\hfill \\ \hfill {?}_{16}^{(k:k+1)}& =& -{?}_{12}^{\left(k\right)}X{P}_R^{J^{k}f}\hfill \\ \hfill {?}_{21}^{(k:k+1)}& =& {?}_{21}^{(k+1)}X{?}_{21}^{\left(k\right)}\hfill \\ \hfill {?}_{22}^{(k:k+1)}& =& {?}_{22}^{(k+1)}-{?}_{21}^{(k+1)}X{?}_{12}^{(k+1)}\hfill \\ \hfill {?}_{23}^{(k:k+1)}& =& {?}_{21}^{(k+1)}XY+{?}_{23}^{(k+1)}\hfill \\ \hfill {?}_{24}^{(k:k+1)}& =& {?}_{21}^{(k+1)}X{?}_{24}^{\left(k\right)}\hfill \\ \hfill {?}_{25}^{(k:k+1)}& =& {?}_{24}^{(k+1)}-{?}_{21}^{(k+1)}X{?}_{14}^{(k+1)}\hfill \\ \hfill {?}_{26}^{(k:k+1)}& =& {?}_{21}^{(k+1)}X{P}_R^{J^{k}f}\hfill \end{array}$$ uid38

The two-handle equations of the resultant assembly express the spatial velocity vectors of the terminal handles of the assembly in terms of the spatial constraint impulses on the same handles, as well as the newly added modal generalized speeds of each constituent flexible body, and the newly introduced $dof$ at the connecting joint. These are the equations which address the dynamics of the assembly when both types of transitions occur simultaneously. In other words, they are applicable when new flexible modes are added to the flexible constituent subassemblies and new degrees of freedom are released at the connecting joint. If there is no change in the joint free-motion map, the spatial partial velocity vector associated with $u^{(k/k+1)f}$ does not appear in the handle equations of the resulting assembly.

5. Hierarchic assembly-disassembly

The DCA is implemented in two main passes: assembly and disassembly [13], [16]. As mentioned previously, two consecutive bodies can be combined together to recursively form the handle equations of the resulting assembly. As such, the assembly process starts at the individual sub-domain level of the binary tree to combine the adjacent bodies and form the equations of motion of the resulting assembly. This process is recursively implemented as that of the binary tree to find the impulse-momentum equations of the new assemblies. In this process, the spatial velocity vector (after transition) and impulsive load of the handles at the common joint of the consecutive bodies are eliminated. The handle equations of the resulting assembly are expressed in terms of the constraint impulses and spatial velocities of the terminal handles, as well as the newly introduce modal generalized speeds and generalized speeds associated with the newly added degrees of freedom at the connecting joints. This process stops at the top level of the binary tree in which the impulse-momentum equations of the entire system are expressed by the following two-handle equations

$$\begin{array}{ccc}\hfill v_1^{1f}& =& {?}_{11}^{(1:n)}{?}_{t_c}^{t_f}{F}_{1c}^{1}dt+{?}_{12}^{(1:n)}{?}_{t_c}^{t_f}{F}_{2c}^{n}dt+{?}_{13}^{(1:n)}\hfill \\ & +& {?}_{14}^{(1:n)}{\stackrel{�}{q}}^{(1:n)f_i}+{?}_{15}^{(1:n)}{u}^{(1:n)f}\hfill \\ \hfill v_2^{nf}& =& {?}_{21}^{(1:n)}{?}_{t_c}^{t_f}{F}_{1c}^{1}dt+{?}_{22}^{(1:n)}{?}_{t_c}^{t_f}{F}_{2c}^{n}dt+{?}_{23}^{(1:n)}\hfill \\ & +& {?}_{24}^{(1:n)}{\stackrel{�}{q}}^{(1:n)f_i}+{?}_{25}^{(1:n)}{u}^{(1:n)f}\hfill \end{array}$$ uid39

Note that through the partial velocity vectors ${?}_{ij}^{(1:n)},(i=1,2\text{and}j=4,5)$, these equations are linear in terms of the newly added generalized modal speeds as well as the generalized speeds associated with the released $dof$ at the joints of the system.

The two-handle equations for the assembly at the primary node is solvable by imposing appropriate boundary conditions. Solving for the unknowns of the terminal handles initiates the disassembly process [17], [18]. In this process, the known quantities of the terminal handles of each assembly are used to solve for the spatial velocities and the impulsive loads at the common joint of the constituent subassemblies using the handle equations of each individual subassembly. This process is repeated in a hierarchic disassembly of the binary tree where the known boundary conditions are used to solve the impulse-momentum equations of the subassemblies, until the spatial velocities of the fine model and impulses on all bodies in the system are determined as a known linear function of the newly introduced generalized speeds of the fine model.

6. Conclusion

The method presented in this chapter is able to efficiently simulate discontinuous changes in the model definitions for articulated multi-flexible-body systems. The impulse-momentum equations govern the dynamics of the transitions when the number of deformable modes changes and the joints in the system are locked or released. The method is implemented in a divide-and-conquer scheme which provides linear and logarithmic complexity when implemented in serial and parallel, respectively. Moreover, the transition from a coarse-scale to a fine-scale model is treated as an optimization problem to arrive at a finite number of solutions or even a unique one. The divide-and-conquer algorithm is able to efficiently produce equations to express the generalized speeds of the system after the transition to the finer models in terms of the newly added generalized speeds. This allows the reduction in computational expenses associated with forming and solving the optimization problem.

Acknowledgment

Support for this work received under National Science Foundation through award No. 0757936 is gratefully acknowledged.