When multiple bodies interact with their neighbours according to a graph.

Introduction

Recently I wrote a code for solving a certain type of multibody problems - just for fun. One instance of such problem is the truss structure, which can be viewed as a set of points connected by bars that can only deform axially. A graph specifies the connectivity of the points - only the points that are connected interact with each other. Another instance of the problem is certain particle system that appear in material science. The particles only interact with their neighbours. Now the problem is essentially determining the equilibrium geometrical configuration of the multi-body system when it is subject to an external excitation. Such equilibrium can be determined via a certain energy principle, if there is any.

Mathematically, the problem is stated as follows. Suppose a multibody system is composed of $N$ points $\{ {\mathbf{x}}_i\}$ with $M$ point pairs $\{ {\mathbf{e}}_i=({\mathbf{x}}_{i1},{\mathbf{x}}_{i2})\}$. The energy ${\mathcal{U}}$ of the multibody system is defined as, $\displaystyle {\mathcal{U}}({\mathbf{u}}) = \sum_{i=1}^M U({\mathbf{x}}_{i1},{\mathbf{x}}_{i2}) \equiv \sum_{i=1}^M U_i$ where $U_i$ is the energy of the $i$th point pair and ${\mathbf{u}}$ is a vector containing all the states of the points. Next, assume the external excitation is conservative, i.e. there exists a scalar function ${\mathcal{V}}({\mathbf{u}})$ such that $-{\frac{\partial {\mathcal{V}}}{\partial {\mathbf{u}}}}$ is the external force acting on the points. The system achieves equilibrium at its stationary point, i.e. when the gradient of total energy $\Pi={\mathcal{U}}+{\mathcal{V}}$ becomes zero, $\displaystyle {\mathbf{g}}({\mathbf{u}}^*) = {\frac{\partial \Pi}{\partial {\mathbf{u}}}} = 0$

Solution schemes

Naive Newton scheme

In simple cases, ${\mathbf{u}}^*$ is a minimum point near the initial state. One can employ any standard optimization algorithms to solve the problem, i.e. to minimize $\Pi$ w.r.t. ${\mathbf{u}}$. Particularly, given the gradient ${\mathbf{g}}$ and the Jacobian ${\mathbf{J}}$ of the energy, $\displaystyle {\mathbf{J}}({\mathbf{u}}^*) = \frac{\partial^2 \Pi}{\partial {\mathbf{u}}^2}$ the minimum is found via the Newton-Raphson (NR) procedure, $\displaystyle \begin{aligned} {\mathbf{J}}_i{\Delta\mathbf{u}}_i &= -{\mathbf{g}}_i \\ {\mathbf{u}}_{i+1} &= {\mathbf{u}}_i + {\Delta\mathbf{u}}_i \end{aligned}$ where ${\mathbf{u}}_0$ is the initial guess of the solution.

By applying the NR procedure to a series of external forces, one obtains the equilibrium path of the system, or the displacement-load curve in structural analysis.

Nonlinear schemes

There are exceptions that cannot be solved using the naive Newton scheme. Such exceptions occur typically when (1) the initial state is far from the solution, or/and (2) the equilibrium path contains limit points. There are two types of limit points. At a load limit point, a nonzero displacement increment results in a zero load increment. Or put in a different way, increasing (or decreasing) the load would result in the increase/decrease of the displacement on the two sides of the limit point. The case of displacement limit point is the reverse.

People have developed various solution schemes to tackle these issues. A review of these schemes is provided in Leon2011. (Also check its Figure 1 for an illustration of the displacement and load limit points.)

To introduce these schemes, in the following we consider only the case where the external excitation is linear, $\displaystyle -{\frac{\partial {\mathcal{V}}}{\partial {\mathbf{u}}}} = \lambda{\mathbf{f}}$ where ${\mathbf{f}}$ is a constant vector and the value of $\lambda$ determines the magnitude of the external force. In such case, the total energy of the system is written as, $\displaystyle \Pi({\mathbf{u}},\lambda) = {\mathcal{U}}({\mathbf{u}}) - \lambda {\mathbf{u}}^T{\mathbf{f}}$

N+1 dimensional space formulation

Assume that the NR iteration starts at the previous equilibrium point. Treating $\lambda$ as an unknown, in the $i$th NR iteration, the linear system is augmented as, $\displaystyle \begin{bmatrix} {\mathbf{J}}_i & -{\mathbf{f}}\\ {\mathbf{c}}_i^T & b_i \end{bmatrix} \begin{bmatrix} {\Delta\mathbf{u}}_i \\ {\Delta\lambda}_i \end{bmatrix} = \begin{bmatrix} -{\mathbf{g}}_i \\ H_i \end{bmatrix}$ where the last row and the coefficients ${\mathbf{c}}_i$, $b_i$ and $H_i$ are introduced to constrain $\lambda$ and close the linear system. Different choices of these coefficients give rise to different solution schemes.

The solution of the augmented system is expressed as, due to Batoz1979, $\displaystyle \begin{aligned} {\Delta\lambda}_i &= \frac{H_i - {\mathbf{c}}_i^T{\Delta\mathbf{u}}_{gi}}{b + {\mathbf{c}}_i^T{\Delta\mathbf{u}}_{fi}} \\ {\Delta\mathbf{u}}_i &= {\Delta\lambda}_i{\Delta\mathbf{u}}_{fi} + {\Delta\mathbf{u}}_{gi} \end{aligned}$ where ${\mathbf{J}}_i{\Delta\mathbf{u}}_{gi}=-{\mathbf{g}}_i$ and ${\mathbf{J}}_i{\Delta\mathbf{u}}_{fi}={\mathbf{f}}$.

Simple incremental schemes

A naive choice of the constraint for $\lambda$ is that $\displaystyle \begin{aligned} {\mathbf{c}}_i=0,\quad b_i=1,\quad H_1={\bar\lambda}\\ H_i=0,\quad i\geq 2 \end{aligned}$ This corresponds to the incremental loading (IL) scheme, which prescribes the load increment ${\bar\lambda}$ at the first iteration. As a result, the load parameter is fixed for the rest of iterations. Near the load limit point, the IL scheme would either fail or result in the snap-through behavior - the solution “jumps” to an equilibrium point far away from the previous solution and results in a non-smooth equilibrium path. Such behavior cannot be avoided unless the direction of the load increment is reversed at the load limit point. Similarly there is incremental displacement scheme which prescribes the displacement increment and suffers from the displacement limit point.

Arc length schemes

The family of arc length (AL) schemes are developed to overcome the limitations of the incremental schemes at limit points. These schemes have been widely used in engineering applications, esp. commercial computater-aided engineering software. The key idea is to introduce an arc length constraint that simultaneously controls the increments in the displacement and load, $\displaystyle {\Delta\mathbf{u}}^T{\Delta\mathbf{u}}+ \mu{\Delta\lambda}^2 = {\Delta s}^2$ Different values of $\mu$ result in the variants of the schemes: $\mu=1$ for spherical AL method, $\mu=0$ for cylindrical AL method, and the most general, $\mu>0,\mu\neq 1$ for elliptical AL method.

For the $N+1$ formulation, the AL constraint is rewritten as, $\displaystyle {\Delta\mathbf{u}}_1^T{\Delta\mathbf{u}}_i + \mu{\Delta\lambda}_1{\Delta\lambda}_i = {\Delta s}_i^2$ where ${\Delta s}_1={\Delta s}$ the prescribed AL and ${\Delta s}_i=0,i\geq 2$. The coefficients ${\mathbf{c}}_i$, $b_i$ and $H_i$ are then determined accordingly. The final expression for load increments is the following, $\displaystyle {\Delta\lambda}_i = \left\{ \begin{array}{ll} \pm\frac{ {\Delta s}}{ \sqrt{ {\Delta\mathbf{u}}_{f1}^T{\Delta\mathbf{u}}_{f1}+\mu}}, &\quad i=1 \\ -\frac{ {\Delta\mathbf{u}}_{f1}^T{\Delta\mathbf{u}}_{gi}}{ {\Delta\mathbf{u}}_{f1}^T{\Delta\mathbf{u}}_{fi}+\mu{\Delta\lambda}_1}, &\quad i\geq 2 \end{array} \right.$

A noteworthy issue of the above scheme is that the sign of ${\Delta\lambda}_1$ is unspecified. One way to determine the sign is to pick the new increment that is “closest to” the previous solution, as discussed in this document. Another way to overcome the issue of direction determination is the linearization of the AL scheme, which replaces the AL constraint with a normal plane constraint.

Generalized displacement control scheme

The Yang1990 scheme is another scheme that bears a similar geometrical interpretation like the linearized AL scheme. Yet, initially, its derivation was to ensure the numerical stability of the $N+1$ dimensional space solution. A generalized stiffness parameter GSP is proposed, $\displaystyle GSP = \frac{ ({\Delta\mathbf{u}}_f^0)^T{\Delta\mathbf{u}}_f^0}{ ({\Delta\mathbf{u}}_{f1}^0)^T{\Delta\mathbf{u}}_{f1}}$ where ${\Delta\mathbf{u}}_f^0$ is the first iteration from the first increment and ${\Delta\mathbf{u}}_{f1}^0$ is the first iteration from the previous increment. At the first increment, ${\Delta\mathbf{u}}_f^0={\Delta\mathbf{u}}_{f1}^0$. Initially $GSP=1$ and GSP increases/decreases as the structure stiffens/softens. A negative value of GSP means that the equilibrium path passed the load limit point and the load direction should be reversed. Refer to Figure 2 in [Yang1990] for a geometrical illustration of the effect of GSP.

The introduction of GSP leads to the following definition of the load increment, $\displaystyle {\Delta\lambda}_i = \left\{ \begin{array}{ll} sgn{\Delta s}\sqrt{ |GSP|}, &\quad i=1 \\ -\frac{ ({\Delta\mathbf{u}}_{f1}^0)^T{\Delta\mathbf{u}}_{gi}}{ ({\Delta\mathbf{u}}_{f1}^0)^T{\Delta\mathbf{u}}_{fi}}, &\quad i\geq 2 \end{array} \right.$ where $sgn=1.0$ and its sign changes whenever GSP becomes negative - and thus the load direction is automatically reversed at load limit point.

Application to geometrically nonlinear truss

Formulation

For the geometrically nonlinear truss, the energy between one pair of points is, $\displaystyle E_i = \frac{1}{2} \frac{k}{||{\mathbf{l}}||} (||{\mathbf{l}}'|| - ||{\mathbf{l}}||)^2$ where ${\mathbf{l}}={\mathbf{x}}_{i2}-{\mathbf{x}}_{i1}$, ${\mathbf{l}}'={\mathbf{x}}_{i2}'-{\mathbf{x}}_{i1}'$, and superscript ’ means new states of the points. The gradient and the Jacobian of $E_i$ only depend on ${\mathbf{u}}_i^T=[{\mathbf{x}}_{i1}^T,{\mathbf{x}}_{i2}^T]$, $\displaystyle \begin{aligned} {\frac{\partial E_i}{\partial {\mathbf{u}}_i}} &= k_1 {\mathbf{L}}\\ \frac{ \partial^2 E_i}{ \partial {\mathbf{u}}_i^2} &= k_1 \begin{bmatrix} I & -I \\ -I & I \end{bmatrix} + k_2 {\mathbf{L}}{\mathbf{L}}^T \end{aligned}$ where, $\displaystyle \begin{aligned} {\mathbf{L}}^T &= [-({\mathbf{l}}')^T, ({\mathbf{l}}')^T] \\ k_1 &= \frac{k}{||{\mathbf{l}}||} \left(1 - \frac{||{\mathbf{l}}||}{||{\mathbf{l}}'||}\right) \\ k_2 &= \frac{k}{||{\mathbf{l}}'||^3} \end{aligned}$ The gradient and Jacobian of the whole system are assembled from those of the point pairs.

Case study

The star dome 24-bar truss is a typical problem that contains various nonlinear behaviors. The description of the problem can be found here. The star dome is loaded at the center and experiences several load and displacement limit points during the loading process. The initial configuration (red dot-dashed line) and several typical deformation are shown in Figure 1. The displacement-load curve of the center node is presented in Figure 2.