Extended system equations
In order to solve the equilibrium equations a equation system with n+1 DOF’s will be formulated. The extra degree of freedom arises from the external load - scale parameter, also denoted as load-proportionality-factor \(\lambda\), which will be appended to the system vector \(\boldsymbol{U}\). The transformation will be shown for both the equilibrium and tangent stiffness. This extension is performed by appending the LPF factor \(\lambda\) to the system displacement column vector \(\boldsymbol{U}\). This new combined system column vector is called \(\boldsymbol{V}\).
(1)\[\begin{split}\boldsymbol{V} = \left[\begin{array}
~ \\
\boldsymbol{U} \\
~ \\ \hdashline
{\lambda}
\end{array}\right]\end{split}\]
Extended system equilibrium
The extended system equilibrium is expanded by one extra equation, which is also referred to as control equation \({g}_{C}(\boldsymbol{V})\).
(2)\[\begin{split}-\boldsymbol{g}_{extended} = \left[\begin{array}
~ \\
-\boldsymbol{g}(\boldsymbol{V}) \\
~ \\ \hdashline
-{g}_{C}(\boldsymbol{V})
\end{array}\right]\end{split}\]
This extra control equation is formulated as a linear constraint to limit the control component j.
(3)\[-{g}_{C} = {V}_j - {V}_{j,max}\]
With this definition the modified system equilibrium becomes
(4)\[\begin{split}-&\boldsymbol{g}_{extended}(\boldsymbol{V}) &= \left[\begin{array}
~ \\
-\boldsymbol{g}(\boldsymbol{V}) \\
~ \\ \hdashline
-{g}_{C}(\boldsymbol{V})
\end{array}\right] = \left[\begin{array}
~ \\
\boldsymbol{r}(\boldsymbol{U}) - \lambda~\boldsymbol{f}_0 \\
~ \\ \hdashline
{V}_j - {V}_{j,max}
\end{array}\right] = \boldsymbol{0} \\
||(-)&\boldsymbol{g}_{extended}(\boldsymbol{V})|| &\le \varepsilon_{tol}\end{split}\]
Extended system tangent stiffness matrix
To illustrate the components of the partial derivatives of the control equation we formulate the first order derivative of the control equation for small changes in \(\boldsymbol{dU}\) and \(d\lambda\) at a given state \((\boldsymbol{U},\lambda)\).
(5)\[\begin{split}-\boldsymbol{g}_{C}(\boldsymbol{V}+\boldsymbol{dV}) &= &-{g}_{C}(\boldsymbol{V}) &- \quad {dg}_{C}(\boldsymbol{V}) &= \\
&= {V}_j &- {V}_{j,max} &- \frac{\partial ({V}_j - {V}_{j,max})}{\partial \boldsymbol{V}}~\boldsymbol{dV} &= 0\end{split}\]
and therefore
(6)\[\frac{\partial {V}_j}{\partial \boldsymbol{V}}~\boldsymbol{dV} = {g}_{C}(\boldsymbol{V})\]
with
(7)\[\begin{split}{dV}_j &= \frac{\partial {V}_j}{\partial \boldsymbol{U}} ~ &\boldsymbol{dU} &+ \frac{\partial {V}_j}{\partial \lambda} ~ &d \lambda \\
{dV}_j &= \boldsymbol{q}_{c,\boldsymbol{U}} ~ &\boldsymbol{dU} &+ {q}_{c,\lambda} ~ &d \lambda\end{split}\]
The modified tangent stiffness \(\boldsymbol{K}_{T,extended}\) is now calculated w.r.t. the modified system vector \(\boldsymbol{V}\).
(8)\[\begin{split}\boldsymbol{K}_{T,extended} = \left[\begin{array}{ccc:c}
~ & ~ & ~ & ~ \\
~ & \frac{\partial \boldsymbol{r}}{\partial \boldsymbol{U}} & ~ & -\boldsymbol{f}_0 \\
~ & ~ & ~ & ~ \\ \hdashline
~ & \frac{\partial {g}_{C}}{\partial \boldsymbol{U}} & ~ & \frac{\partial {g}_{C}}{\partial \lambda}
\end{array}\right]\end{split}\]
(9)\[\begin{split}\boldsymbol{K}_{T,extended} = \left[\begin{array}{ccc:c}
~ & ~ & ~ & ~ \\
~ & \boldsymbol{K}_{T} & ~ & -\boldsymbol{f}_0 \\
~ & ~ & ~ & ~ \\ \hdashline
~ & \boldsymbol{q}_{C,\boldsymbol{U}} & ~ & {q}_{C,\lambda}
\end{array}\right]\end{split}\]
with the partial derivatives of the control equation:
(10)\[\begin{split}\boldsymbol{q}_{C,\boldsymbol{U}} &= \frac{\partial {g}_{C}(\boldsymbol{U},\lambda)}{\partial \boldsymbol{U}} \\
\boldsymbol{q}_{C,\lambda} &= \frac{\partial {g}_{C}(\boldsymbol{U},\lambda)}{\partial \lambda}\end{split}\]
The combined partial derivate of the control equation may be expressed as:
(11)\[\boldsymbol{q}_{C} = \begin{bmatrix}
\boldsymbol{q}_{C,\boldsymbol{U}} & {q}_{C,\lambda}
\end{bmatrix}\]
and is a system row vector with a length of nDOF+1 entries filled with zeros, except for the j-th entry (control component), which is one. Let’s assume we have a system of 6 DOF and the control component is identified to j=4 then \(\boldsymbol{q}_{C}\) looks like
(12)\[\boldsymbol{q}_{C} = \begin{bmatrix}
0 & 0 & 0 & 1 & 0 & 0 & 0
\end{bmatrix}\]
Summary
The final equation system for the Path-Following Algorithm may now be formulated as
(13)\[\boldsymbol{K}_{T,extended} ~ \boldsymbol{dV} = \boldsymbol{g}_{extended}(\boldsymbol{V})\]
and in detail
(14)\[\begin{split}\begin{bmatrix}
\boldsymbol{K}_{T} & -\boldsymbol{f}_0 \\
\boldsymbol{q}_{c,\boldsymbol{U}} & \boldsymbol{q}_{c,\lambda}
\end{bmatrix} \begin{bmatrix}
\boldsymbol{dU} \\
d{\lambda}
\end{bmatrix} = \begin{bmatrix}
\boldsymbol{g} \\
{g}_{C}
\end{bmatrix}\end{split}\]
The control component \(j\) is defined as the signed index of the biggest component of the incremental system vector \(\boldsymbol{dV}\). This component remains fixed during all newton-iterations inside an increment. To initialize the control component \(j\) the linear equation system is solved with an incremental load proportionality factor \(d\lambda\).
(15)\[\begin{split}\boldsymbol{K}_{T} ~ \boldsymbol{dU} &= d \lambda \boldsymbol{f}_0 \\
\rightarrow \boldsymbol{dU} &= \dots\end{split}\]
(16)\[j = \text{index} [\max (|\boldsymbol{dU}|)] \cdot \text{sign} [\max (|\boldsymbol{dU}|)]\]