Solvers and eigen value decomposition

To solve the structural forward dynamics of the system \eqref{eq:C5:femmodel}, the toolkit uses an implicit Newmark-\(\beta\) solver \cite{Newmark1959Jul}, which is briefly outlined in Appendix \ref{app:C5:newmark}. \(\renewcommand{\elastic}{_\textrm{e}}\) \(\renewcommand{\grav}{_\textrm{g}}\) \(\renewcommand{\p}{\partial}\) \(\renewcommand{\R}{\mathbb{R}}\) \(\renewcommand{\vec}[1]{\mathbf{#1}}\) \(\renewcommand{\mat}[1]{\mathbf{#1}}\) \(\renewcommand{\FB}{\mathbf{F}}\) \(\renewcommand{\fB}{\mathbf{f}}\) \(\renewcommand{\Hm}{\mathcal{H}}\) \(\renewcommand{\xB}{\mathbf{x}}\) \(\renewcommand{\AB}{\mathbf{A}}\) \(\renewcommand{\x}{\mathbf{x}}\) \(\renewcommand{\dxB}{\dot{\mathbf{x}}}\) \(\renewcommand{\dt}{\delta t}\) \(\renewcommand{\ddxB}{\ddot{\mathbf{x}}}\)

Explaination on Newmark-\(\beta\) solver

The Newmark-\(\beta\) method is an implicit numerical integration scheme extensively used to solve high-dimensional structural dynamic problem ¹². We briefly explain the algorithm implemented in the function Fem.simulate. First, let us subdivide the time domain such that \((0,...,T)\) with uniform timesteps \(\dt = t_{i+1} - t_i\). Then, given the initial conditions, we wish to compute the state evolution \(\xB(t_i)\) and \(\dot{\x}(t_i)\). For conciseness, let us write the discrete states of the FEM model as \(\xB(t_i) = \xB^{(i)}\). Through the extended mean value theorem, we can formulate the general Newmark-\(\beta\) scheme as

\[ \begin{align} {\dxB}^{(i+1)} & = {\dxB}^{(i)} + \dt \left[(1-\beta_1)\ddxB^{(i)} + \beta_1 \ddxB^{(i+1)} \right], \\[0.2em] % {\xB}^{(i+1)} & = {\vec{x}}^{(i)} + {\dt}\left[ \dxB^{(i)} +\dt(\tfrac{1}{2}-\beta_2) \ddxB^{(i)} + \dt \beta_2 \ddxB^{(i+1)} \right], \end{align} \]

where \(\beta_1,\beta_2 \ge \frac{1}{2}\). Now, in the expressions above only the forward-time acceleration \(\ddxB^{(i+1)}\) is the unknown partial solution, hence we conveniently write \(\vec{w} := \ddot{\xB}^{(i+1)}\). Substitution into the flow \eqref{eq:C5:femmodel}, we find:

\[ \mat{r}(\vec{w}): = \mat{M} \vec{w} + \nabla_{\xB\,}\mathcal{U}(\vec{w}) + \mat{R}\dxB^{(i+1)}(\vec{w}) - \vec{G}\vec{u}^{(i+1)}, \label{app:C3:residual_newmark} \]

where \(\Hm\) is the Hamiltonian. Following, the residual vector \eqref{app:C3:residual_newmark} forms an optimization problem in the form \(\text{argmin}_{\vec{w}} \lVert \vec{r}(\vec{w}) \rVert_2\) for unknown accelerations \(\vec{w}\). This implicit relation can be solved numerically using a recursive Newton Raphson method. Given the \(n\)-th iteration, the recursive solver reads

\[ \vec{w}^{(n+1)} = \vec{w}^{(n)} - \alpha_+ \left[ \mat{A}(\vec{w}^{(n)}) \right]^{-1}\! \vec{r}(\vec{w}^{(n)}), \label{eq:C3:newmark_newtonsolver} \]

where \(\AB := \left[\mat{M} + \beta_1 \dt \mat{R} + \beta_2 \dt^2 \mat{K}_{T}\right]\) is the hessian matrix, and \(0< \alpha_+ \le 1\) an update coefficient. The matrix \(\mat{K}_T\) denotes the tangent stiffness related to the local gradient of the elasticity force, given by \(\mat{K}_T := \nabla_{\xB} \fB\elastic\).

Implicit solvers offer improved stability compared to explicit methods, such as the Runge-Kutta solver (\texttt{ode45}), particularly when larger time steps are employed. However, the cost of larger time steps is a decreased numerical precision. Alternatively, for quasi-static problems when \(\ddxB = \dxB = \vec{0}_n\), we aim to seek the solutions to the static force equilibrium \(\vec{r}(\x) = \vec{0}_n\) where \(\vec{r} := -\fB_{\textrm{mat}} + \fB\grav + \fB_\textrm{u} + \fB_{\Omega_{\textrm{env}}}\) is the force residual vector. The nonlinear equality for nodal displacements \(\x\) is solved using an iterative Newton-Raphson solver. To call these solvers, dynamic simulations are executed with \code{fem.simulate()} and quasi-static simulations with fem.solve(). Upon completion of a simulation, all displacements, velocities, forces, and stress information are stored in the \code{fem.Log} data structure. This log file can be accessed for data analysis or during simulation to facilitate state feedback control.

Example: basic pull test

% generate mesh
msh = Mesh(sRectangle(20), 'Quads', [25,10]);
msh = msh.generate();

% generate fem model and boundary cond.
fem = Fem(msh, 'TimeStep', 1/100);
fem = fem.addSupport('left', [1, 1]);
fem = fem.addDisplace('right', [100, 0]);
fem = fem.addMaterial(NeoHookean(1, 0.33));

% solve :)
fem = fem.solve();

Alternatively, we can explore nonlinear modal analysis at any quasi-static equilibrium configuration \(\x^* \in \mathcal{X}\) of the system. Let \(\mat{K}_T:= \left[\frac{\p \fB\elastic}{\p x_1} \;...\; \frac{\p \fB\elastic}{\p x_2} \right]\) be the Jacobian matrix of the (nonlinear) elastic potential forces, also referred to as the tangent stiffness. Then, the local eigenvalue problem for the linearized FEM model around the point \(\x^*\) is given by

\[ \Big[\mat{K}_T(\x^*)- \lambda_i \mat{M} \Big] \vec{\theta}_i = \vec{0}_n, \]

where \(\lambda_i\) is a real scalar eigenvalue and \(\vec{\theta}_i\) is its corresponding eigenmode. The dynamic analysis is implemented in Sorotoki using \code{fem = fem.analysis(x)}, which stores the necessary data in \code{fem.Log}. It is important to note that, unlike linear finite element models, the set of eigenmodes \({\vec{\theta}_i}\) obtained from the eigenvalue decomposition in \eqref{eq:C5:eigenmode_decomposition} is highly dependent on the linearization point \(\x^*\) and may thus not be unique for all \(\x^* \in \mathcal{X}\).

Nathan M. Newmark. A Method of Computation for Structural Dynamics. Journal of the Engineering Mechanics Division, 85(3):67–94, 1959. doi:10.1061/JMCEA3.0000098. ↩
Gerhard A. Holzapfel. Nonlinear Solid Mechanics: A Continuum Approach for Engineering Science. Volume 37. Kluwer Academic Publishers, Heidelberg, Germany, July 2002. doi:10.1023/A:1020843529530. ↩