The topics cross-cuts several research activities, and it concerns methodological aspects of the numerical strategies employed in structural dynamics problems.
Recently, together with my colleague Dr. Franco Milicchio we achieved a significant milestone by developing in C++ a quite general-purpose for path-following periodic solutions of generically non-smooth dynamics systems. We published the fully open-source code on Zenodo:
Improving the monodromy matrix computation in pathfollowing schemes for nonsmooth dynamics. International Journal of Non-Linear Mechanics 2023, DOI: 10.1016/j.ijnonlinmec.2023.104455 CONTRIBUTORS: Formica, G.; Milicchio, F.; Lacarbonara, W.
We enhanced the pseudo-arclength path-following scheme proposed in 2022 by focusing on the numerical computation of the monodromy matrix, one of the key points of the numerical strategy. Such matrix, by serving as iteration matrix, governs indeed both the accuracy and robustness of the whole strategy.
In nonsmooth dynamic systems, it can only be computed via numerical differentiation (eg, central difference scheme), which leads to subtractive cancellation errors.
We proposed a new technique by performing a Lanczos differentiation by integration of the vector field associated with the ODE problems, directly within the time-step integration. In this way, we are able to bound the errors induced by the standard central difference scheme.
A Runge-Kutta single-pass method is then implemented to construct the monodromy matrix
We conducted a numerical campaign of tests with several degrees of freedom, by accounting for both hardening-type and softening-type nonlinearities.
We also compared the results with the analytical Jacobian, whenever this is available.
The numerical tests showed that our technique allows the algorithm to achieve a higher level of accuracy when compared with standard finite differentiation employed after time integrations of the perturbed periodic solutions.
In a 20-dimensional nonlinear hardening system, we were able to reduce the mean absolute error by four orders of magnitude.
Additionally, for hysteretic systems we proved that our approach is capable of providing better robustness to the whole numerical strategy, and of improving its convergence properties. In a 30-dimensional hysteretic system, our single-pass method was able to attain convergence employing a step size of the solution path four times larger than the standard method.
Finally, the CPU time per step size is about 3.5 lower than the CPU time required by the previously proposed algorithm.
The resulting C++ code appears as a 2nd version at our
Zenodo web repo
A Krylov accelerated Newton–Raphson scheme for efficient pseudo-arclength pathfollowing. International Journal of Non-Linear Mechanics 2022, DOI: 10.1016/j.ijnonlinmec.2022.104116 CONTRIBUTORS: Formica, G.; Milicchio, F.; Lacarbonara, W.
We tackled here the computational efficiency of a pseudo-arclength path-following scheme tailored for general multi-degree-of-freedom (multi-dof) nonlinear dynamical systems. The path-following approach is based on the numerical computation of the Poincaré map and its Jacobian in order to treat nonautonomous systems with discontinuous vector fields.
The scheme is applied to obtain frequency response curves of multi-dof hysteretic systems with a state vector size up to 120, as well as various reduced-order models of single and multiple cantilever beams on a shuttle mass.
We proposed an enhanced method for drastically increasing the speed of convergence in the modified Newton–Raphson scheme which is supported by a Krylov sub-space iteration. We also adopted the LU decomposition of a frozen Jacobian matrix, which, upon convergence, becomes the monodromy matrix.
Several numerical tests performed on mechanical systems with material or geometric nonlinearities corroborate the efficiency of the numerical strategy.
The implemented C++ code runs with low memory usage requirements, here 6–10 MB. It was shown that our approach saves over 90% of the CPU with respect to the standard NR scheme.
Although not parallelized, the software is efficient on every platform, and easily portable on every operating system. As a proof of its portability, the algorithm was successfully compiled within an ARM CPU (see incourant project).
The obtained results paves the way towards possible future improve- ments, such as parallelizability of some processes, as well as coding of some computations to be processed in GPU for a generalized speedup of the entire algorithm.
Pathfollowing of high-dimensional hysteretic systems under periodic forcing. Nonlinear Dynamics 2021, DOI: 10.1007/s11071-021-06374-7 CONTRIBUTORS: Formica, G.; Vaiana, N.; Rosati, L.; Lacarbonara, W.
We investigated here some computational aspetcs when the dynamic response and bifurcations of high-dimensional systems endowed with hysteretic restoring forces in all degrees of freedom have to be analyzed.
Two types of hysteresis models are considered, namely the Bouc–Wen model and a differential version of the so-called exponential model of hysteresis.
Both models are generally capable of representing accurately hysteresis phenomena, typical of mechanical systems, devices and materials that exhibit rate-independent hysteretic behavior.
Compared with the former, the latter offers some advantages for an easy use in research and practice: it requires a set of parameters having a clear mechanical meaning; it can be calibrated by means of identification procedures which are simpler than those typically required by the Bouc-Wen models; finally, the exponetial model can be generalized to reproduce asymmetric hysteresis loops.
For both models, the numerical technique tailored for tackling high-dimensional hysteretic systems is here based on an enhanced pathfollowing approach based on the Poincaré map. The numerical strategy uses several computational improvements such as the simplified Newton–Raphson method in which the Jacobian of the map is computed at the initial guess and is not further updated at each iteration and the computation of the columns of the matrix is performed using parallel computing.
The obtained results in terms of frequency-response curves together with their stability show the richness of responses of hysteretic systems in terms of switch of nonlinearity depending on the constitutive parameters, a richer multi-stability and, in general, a significant sensitivity to the constitutive hysteresis parameters.
In particular, a five-dof mass-spring-damper-like system, with each rheological element described by the Bouc– Wen or the exponential model of hysteresis enriched by cubic and quintic nonlinear elastic terms, is inves- tigated and a rich variety of nonlinear responses and bifurcations is found and discussed.
In the Figure, the set of frequency-response curves, together with some force–displacement hysteresis loops, shows the pronounced softening effect and the appearance of multi-stability regions due to the fold bifurcations for some value of the nondimensional base acceleration amplitude. For other values, there is a coalescence of the twofold bifurcations into a co-dimension-two bifurcation point past which the fold points disappear and the response becomes stable in the whole range of primary resonance frequencies.
Computational efficiency and accuracy of sequential nonlinear cyclic analysis of carbon nanotube nanocomposites. Advances in Engineering Software 2018, DOI: 10.1016/j.advengsoft.2018.08.013 CONTRIBUTORS: Formica, G.; Milicchio, F.; Lacarbonara, W.
Designing new multiphase, multiscale nanostructured materials via a simulation-driven approach is extremely challenging due to the high computational cost and time. Especially for capturing nonlinear and hysteretic behavior, several optimization cycles are required, keeping both efficiency and accuracy.
For this puropose, the computational approach resorting to the nonlinear 3D finite element implementation uses a variant of the Newton-Raphson method within an enhanced time integration scheme.
Such an approach is inspired to the work of my mentor Prof. Casciaro
(DOI: 10.1007/BF02149027).
In particular, it works with the elastic tangent matrix as iteration matrix, so as to avoid its iterative update which higly expensive from a computational standpoint. The single iteration matrix assembled just once per both the single cyclic analysis and the overall damping curve.
This is especially rewarding in the context of the employed mechanical model which exhibits hysteresis manifested through a discontinuous change in the stiffness at the reversal points where the loading direction is reversed.
Key implementation aspects – such as the integration of the nonlinear 3D equations of motion, the numerical accuracy/efficiency as a function of the time step or the mesh size – were discussed in this paper. In particular, efficiency is regarded as performing fast computations especially when the number of cyclic analyses becomes large. By making use of laptop CPU cores, a good speed of computations is achieved not only through parallelization but also employing a caching procedure for the iteration matrix. Preliminary tests show that our approach has a very low sensitivity to both mesh and time step size refinements. An extensive numerical campaign revealed that just with the serial code, the caching mechanism allows one to obtain a good speedup in the damping curve with a quasi-linear behavior. Furthermore, when executing the parallelized code without using cache, we were able to scale well with the number of employed cores in an off-the-shelf mini-desktop computer. This performance turned out to be improved by coupling the effects of our parallel execution and caching mechanism. The findings of this research pave the way towards fast genetic-like optimizations of nanocomposite multilayered structures subject to complex multi-axial loading conditions.
Coupling FEM with parameter continuation for analysis of bifurcations of periodic responses in nonlinear structures. Journal of Computational and Nonlinear Dynamics 2013, DOI: 10.1115/1.4007315 CONTRIBUTORS: Formica, G.; Arena, A.; Lacarbonara, W.; Dankowicz, H.
This was my first step in investigating enhanced approaches for numerical strategies in structural dynamics.
It was taken during my visit to University of Illinois at Urbana-Champaign, that I spent in november 2010 thanks to Prof. Harry Dankowicz.
In the paper, a computational framework was proposed to perform parameter continuation of periodic solutions of nonlinear, distributed-parameter systems represented by partial differential equations with time-dependent coefficients and excitations.
The path-following procedure, encoded in the general-purpose MATLAB-based computational continuation core (COCO), employs only the evaluation of the vector field of an appropriate spatial discretization; for example as formulated through an explicit finite-element discretization or through reliance on a black-box discretization.
An original contribution of this paper was a systematic treatment of the coupling of COCO with COMSOL MULTIPHYSICS, demonstrating the great flexibility afforded by this computational framework. COMSOL MULTIPHYSICS provides embedded discretization algorithms capable of accommodating a great variety of mechanical/physical assumptions and MULTIPHYSICS interactions.
Within this framework, we showed that a concurrent bifurcation analysis may be carried out together with parameter continuation of the corresponding monodromy matrices. As a case study, we considered a nonlinear beam, subject to a harmonic, transverse direct excitation for two different sets of boundary conditions and demonstrate how the proposed approach may be able to generate results for a variety of structural models with great ease. The numerical results include primary-resonance, frequency-response curves together with their stability and two-parameter analysis of multistability regions bounded by the loci of fold bifurcations that occur along the resonance curves.
As the above Figure shows, we investigated the influence of the number of spatial mesh intervals in the finite-element discretization on the evaluation of the nonlinear dynamic response of the beam subject to a primary resonance of the second, third, and fourth symmetric modes. At high excitation amplitudes (in the figures, the lower excitation ampli- tude is 10k0 and the higher excitation amplitude is 50k0), the geo- metric nonlinearities induce significant nonlinear spatial distortions (spatial over-tones) with respect to the linear mode shapes of the excited mode. A refined, finite-element discretization is thus needed to describe correctly the dynamic behavior of the beam at such excitation magnitudes.
We also tested as various boundary conditions, which are straightforward to implement using the appropriate COMSOL commands, lead to different nonlinear contributions to the beam response. The curves at the bottom-right highlight a remarkable effect of changing boundary conditions, where the response switches from hardening to softening as already suggested in previous theoretical and experimental investigations reported in the literature.