Non-intrusive reduced order modelling for the dynamics of geometrically nonlinear flat structures using three-dimensional finite elements
Alessandra Vizzaccaro, Arthur Givois, Pierluigi Longobardi, Yichang Shen, Jean-François Deü, Loïc Salles, Cyril Touzé, Olivier Thomas
NNon-intrusive reduced order modelling for the dynamics of geometricallynonlinear flat structures using three-dimensional finite elements
Alessandra Vizzaccaro · Arthur Givois · PierluigiLongobardi · Yichang Shen · Jean-Fran¸cois De¨u · Lo¨ıcSalles · Cyril Touz´e · Olivier ThomasAbstract
Non-intrusive methods have been used since two decades to derive reduced-order models for geometricallynonlinear structures, with a particular emphasis on the so-called STiffness Evaluation Procedure (STEP), relyingon the static application of prescribed displacements in a finite-element context. We show that a particularly slowconvergence of the modal expansion is observed when applying the method with 3D elements, because of nonlinearcouplings occurring with very high frequency modes involving 3D thickness deformations. Focusing on the case offlat structures, we first show by computing all the modes of the structure that a converged solution can be exhibitedby using either static condensation or normal form theory. We then show that static modal derivatives provide thesame solution with fewer calculations. Finally, we propose a modified STEP, where the prescribed displacements areimposed solely on specific degrees of freedom of the structure, and show that this adjustment also provides efficientlya converged solution.
Keywords
Reduced order modeling · Geometric nonlinearities · Three-dimensional effect · Thickness modes · Modified STiffness Evaluation Procedure · Nonlinear modes · Modal derivatives
Geometrically nonlinear effects appear generally in thin structures such as beams, plates and shells, when the amplitudeof the vibration is of the order of the thickness [26,37]. The von K´arm´an family of models for beams, plates and shellsallows one to derive explicit partial differential equations (PDE) [18,1,37,39,8], showing clearly that a coupling betweenbending and longitudinal motions causes a nonlinear restoring force of polynomial type in the equations of motion. This
A. VizzaccaroVibration University Technology Centre, Imperial College London, SW7 2AZ London, UKE-mail: [email protected]. Givois, O. ThomasArts et Metiers Institute of Technology, LISPEN, HESAM Universit´e, F-59000 Lille, FranceA. Givois, J.-F. De¨uStructural Mechanics and Coupled Systems Laboratory, Conservatoire National des Arts et M´etiers, 2 Rue Cont´e, 75003 Paris, FranceP.Longobardi, L. SallesVibration University Technology Centre, Imperial College London, SW7 2AZ London, UKY. Shen, C. Touz´eIMSIA, ENSTA Paris, CNRS, EDF, CEA, Institut Polytechnique de Paris, 828 Boulevard des Mar`echaux, 91762 Palaiseau Cedex,France a r X i v : . [ c s . C E ] S e p A. Vizzaccaro, A. Givois, P. Longobardi, Y. Shen, J.-F. De¨u, L. Salles, C. Touz´e, O. Thomas geometric nonlinearity is then at the root of complex behaviours, that also need dedicated computational strategies inorder to derive quantitative predictions. On the phenomenological point of view, structural nonlinearities give rise tonumerous nonlinear phenomena that have been analysed in a number of studies: frequency dependence on amplitude[25,20,14], hardening/softening behaviour [46,45], hysteresis and jump phenomena [24,36], mode coupling throughinternal resonances [24,39,9], bifurcations and loss of stability [44,10], chaotic and turbulent vibrations [5,3]. On thecomputational point of view, nonlinear couplings break the invariance property of the linear eigenmodes. Consequently,deriving reduced-order models (ROM) is no longer a straightforward problem and care has to be taken in order to findout a ROM that is capable of describing the dynamics of the whole system without losing accuracy.When the structure under study is discretized with the finite element (FE) method, the problem of deriving accurateROM is more stringent since the user cannot rely on a PDE in order to unfold an ad hoc mathematical method forbuilding the ROM. Moreover, because of the intrinsic nature of the geometrical nonlinearities, all degrees of freedomof the FE model are nonlinearly coupled and substructuring techniques are not suitable, on the contrary to localizednonlinearities occurring frequently in contact and friction problems [17,52]. Consequently, for geometric nonlinearity,the computation of the nonlinear coupling coefficients is as important as finding out a correct reduced basis, andsometimes the two problems are interwoven. Moreover, a number of recent studies highlighted the importance ofusing so-called indirect or non-intrusive methods, where the idea is to use the standard operations that any FE codeis used to perform in order to build the ROM, hence avoiding the need to enter in the code and write new linescomputing the needed quantities. On the other hand, direct methods also exist, for which there is a need to implementthe computations at the level of the element [32,47]. The STEP (STiffness Evaluation Procedure) is a non-intrusivemethod and has been first introduced by Muravyov and Rizzi in 2003 [23]. In its first version as described in [23], itallows computation of the nonlinear coupling coefficients of the discretized problem in the modal basis from a seriesof static computations with prescribed modal displacements. It has then been used in a number of contexts [21,19,22,28,8], and is generally connected to the modal basis. However one has to understand that per se , STEP is just anevaluation technique, a computational non-intrusive method, that can be used with other inputs than those from themodal basis.The STEP, although being largely applied in numerous cases, is known to suffer from a number of problems,making it not as so simple as its formulation could let one think. A first one lies in the amplitude of the prescribeddisplacement one has to impose in order to excite sufficiently the nonlinearity. As shown in [8], there is a clear range ofamplitude for which the method works properly, between a minimal value where nonlinearity is not sufficiently excitedand a maximum value for which other nonlinear effects are appearing. Another problem is connected to the use of themodal basis with a STEP computation, and must be interpreted as a drawback of using the modal basis for nonlinearcomputation, but is not directly linked with the STEP calculation. The problem is that of the slow convergence linkedto the loss of invariance of eigenspaces, and the numerous couplings between low frequency bending modes and highfrequency longitudinal modes, as underlined in a number of papers, see e.g. [22,8] and references therein. Consequently,numerous methods have been proposed in order to overcome this limitation: dual modes [16], POD modes [29], discreteempirical interpolation [41], modal derivatives [12,51], quadratic manifold [13,30], just to name a few.Incidentally, the majority of papers where the STEP has been applied to thin structures, report results obtainedwith structures discretized with beam, plate or shell finite elements, as being mostly interested in slender structures. A few examples with block elements can be found in the literature, see e.g. [27,50], where dual modes have been usedin order to achieve convergence. Using 3D finite elements can be interesting in practice, sometimes mandatory. Inengineering applications, structures are often defined with a 3D geometrical model, for which a 3D FE discretizationis straightforward. In other cases, some particular physical effects (piezoelectricity for instance) are sometimes imple-mented in existing FE codes only with 3D elements. Preliminary studies of the authors reveal a number of difficultieswhen blindly applying the STEP with modal basis to 3D finite elements, with more stringent inaccuracies and prob-lems than those encountered with plate or shell elements. In particular, an unexpected slow convergence was observedwhen using the modal basis, and very high frequency modes appear to be involved.The objective of this paper is to diagnose properly the issues one can encounter when applying the STEP with amodal basis to a structure discretized with 3D finite elements and present methods to overcome the problems. In the on-intrusive ROM for the dynamics of geometrically nonlinear flat structures using 3D FEM 3 course of the paper, we will show that the problem is intrinsically related to the use of the modal basis, and that theSTEP can be used with other inputs than the modal basis, in order to get better results. We restrict our attention tothin structures that are symmetric in the thickness direction, such as straight beams or plates, for which there is abending / longitudinal uncoupling at the linear level that greatly simplifies the understanding of the phenomena andenables to obtain reference results.The paper is organized as follows. Section 2 is dedicated to the framework of the study, the equations of motion anda brief recall of the STEP. Section 3 addresses on test examples the issues of the STEP applied to 3D FE. It is shownthat a nonlinear coupling of bending modes with very high frequency modes involving deformations in the thicknessof the structures occurs, and thus called thickness modes. Those modes are the result of 3D deformation effects thatare not present when using beam or plate models. This unexpected coupling is different from the traditional bending-longitudinal coupling and the numerical examples show that they are of prime importance to achieve a convergedsolution. If one can compute all the coupled modes, two strategies are given to overcome the large dimension of thereduced basis: static condensation and the reduction to a single nonlinear normal modes. Both method shows that whena single master mode drives the dynamics, the ROM can still be composed of a single nonlinear oscillator. However,finding out all the coupled high-frequency modes is generally out of reach for complex structures. In section 4, we thenshow that using a static condensation of a single modal derivative allows to retrieve the same converged result, in amore direct and efficient way. In addition, a modified version of the STEP is proposed, to directly embed the couplingwith thickness modes. It consists in applying the prescribed displacement only on selected degrees of freedom and letthe other free. In section 5, the physical mechanisms of those nonlinear couplings are explained and section 6 presentsnumerical results to validate the proposed numerical methods, able to overcome the convergence issues of the classicalSTEP. x gathers all the degrees of freedom of the model (displace-ments/rotations at each nodes) and is N -dimensional. The equation can be written: M ¨ x + C ˙ x + f ( x ) = f e , (1)where M and C are the N × N dimensional mass and damping matrices, f ( x ) is the internal force vector, f e ( t ) is theexternal force vector and the overdot represents the classical differentiation with respect to time t : ˙ • = d • / d t . Notethat since we are more interested in the computation of the nonlinear restoring force, a linear viscous damping modelhas been used. In the present case of geometrically nonlinear structures, the internal force vector encompasses onlypolynomial terms up to order three and involves the displacement vector x only [19,22,47,11]. For the present study,we restrict our attention to this case, but extensions to more complex cases involving for examples velocity terms can also be handled.It is convenient to split the internal force vector into a linear part and a purely nonlinear part. Assuming that theequilibrium point, the structure at rest, is given by x = , the tangent stiffness matrix K classically writes: K = ∂ f ∂ x (cid:12)(cid:12)(cid:12)(cid:12) x = , (2)so that the nonlinear internal force vector is defined as f nl ( x ) = f ( x ) − Kx , (3) A. Vizzaccaro, A. Givois, P. Longobardi, Y. Shen, J.-F. De¨u, L. Salles, C. Touz´e, O. Thomas and the equations of motion reads: M ¨ x + C ˙ x + Kx + f nl ( x ) = f e , (4)where the geometrically nonlinear part of the problem is concentrated in the nonlinear internal force vector f nl ( x ).The modal basis can be used as a first step in order to make the linear terms diagonal. The eigenmodes are thefamily of couples eigenfrequency-eigenvectors ( ω k , φ k ), k = 1 , . . . , N , solutions of the undamped, free and linearizedEq. (4): ( K − ω k M ) φ k = . (5)Assuming the modal expansion for the displacement vector: x ( t ) = N (cid:88) k =1 φ k X k ( t ) , (6)where X k ( t ) is the modal amplitude, and using Galerkin projection allows one to rewrite the equations of motion inthe modal space, for all k = 1 , . . . , N , as:¨ X k + 2 ζ k ω k ˙ X k + ω k X k + N (cid:88) i =1 N (cid:88) j = i α kij X i X j + N (cid:88) i =1 N (cid:88) j = i N (cid:88) l = j β kijl X i X j X l = Q k , (7)with Q k = φ T k f e /m k , m k = φ T k M φ k , (8)where m k is the modal mass, and where a modal damping (of factor ζ k ) has been assumed uncoupled (an assumptionvalid for small damping, even with non-proportional C matrix [7]). The nonlinear part of this reduced order is writtenwith quadratic and cubic polynomial terms with coefficients α kij and β kijl . It is an exact expansion in the case of 3DFE [47] or beam/plate/shell FE based on von K´arm´an strain/displacement law [19], whereas it is truncated in thecase of geometrically exact theories [38].In Eqs. (7), the linear parameters ω k , φ k and Q k are obtained by the modal analysis of Eq. (5), available in anyfinite element code. The main issue is thus to compute the nonlinear coupling coefficients α kij and β kijl . The STEP(STiffness Evaluation Procedure) when used with the modal basis as first introduced by Muravyov and Rizzi [23], is anon-intrusive (or indirect) method, allowing one to get these coupling coefficients from standard computations availablein any FE code. It relies on imposing prescribed static displacements having the shapes of selected eigenmodes, witha given amplitude. From a clever choice of the modes and the amplitudes, a simple algebra allows one to retrieveall the coefficients from the internal force vector given by the FE code, the key idea being to impose plus/minus thedisplacement with selected combinations of modes. The reader can find detailed explanation on this calculation ina number of papers, including the original one [23], as well as the improvement proposed in [28] using the tangent stiffness matrix.To illustrate the method, we only show the computation of coefficients α kpp and β kppp ; for the general case theinterested reader is referred to [23,28]. The following static displacements are prescribed to the structure: x p = ± λ φ p ⇒ (cid:26) q p = λ,q j = 0 ∀ j (cid:54) = p, (9)where λ refers to an amplitude coefficient of the eigenvector φ p whose value has to be chosen so as to activate thegeometrical nonlinearities. A value of h/
20, where h is the thickness of the plate, is recommended in [8]. Since themodes are orthogonal, imposing a displacement along mode p is equivalent to consider that only the modal coordinate on-intrusive ROM for the dynamics of geometrically nonlinear flat structures using 3D FEM 5 q p is not vanishing, as detailed in the second part of Eq. (9). Since x p is time independent, introducing Eq. (9) intoEqs. (4) and (7) leads to, for all k = 1 , . . . , N : λ α kpp + λ β kppp = φ T k f nl ( λ φ p ) /m k , (10a) λ α kpp − λ β kppp = φ T k f nl ( − λ φ p ) /m k . (10b)Hence the unknown quadratic and cubic coefficients are found easily as α kpp = 12 λ φ k m k ( f nl ( λ φ p ) + f nl ( − λ φ p )) , (11a) β kppp = 12 λ φ k m k ( f nl ( λ φ p ) − f nl ( − λ φ p )) . (11b)Similar algebraic manipulations with more modes involved in the prescribed displacement then allows one to get thefull family of quadratic and cubic coefficients. The non intrusive nature of the method appears clearly: in the FEsoftware, one prescribes the displacement field x p and computes the corresponding external force vector f e = f ( x p )with successive linear and nonlinear static computations. Then, f nl ( x p ) is obtained with Eq. (3).2.2 The case of flat structuresIn this article, we restrict the analysis to flat structures, such as straight beams or plates with a boundary of arbitraryshape. The thickness of the plate can be non-constant, but the geometrical and material distribution must be symmetricwith respect to the middle line / plane of the structure.If a plate theory, with a Kirchhoff-Love kinematics for instance, is applied to this structure, the displacement fieldof any point of the structure is described by the displacement field of the middle plane. There is a membrane / bendingdecoupling in linear elasticity and two families of modes are obtained: bending modes, for which only the transversecomponent of the middle plane displacement are non-zero, and membrane modes, for which the transverse componentof the middle plane displacement field is zero. Analogous properties are valid for a beam theory, with a longitudinal /transverse decoupling on the middle line.In the present case of a 3D structure which has the shape of a beam / plate, it is also possible to split theeigenmodes into two families in the same manner. The first family includes the bending modes, analogous to the onesobtained in the plate theory. Their frequencies are in the lower part of the spectrum, while their deformed shapesare dominated by transverse displacements. They, up to the accuracy of the plate theory, have the same displacementfield in the middle plane / line, with no longitudinal displacement. The second family gathers all the other modes,denoted by non-bending (NB) modes, that appear at higher frequencies in the spectrum. Some of them are analogousto the longitudinal modes of the plate theory, with the same transverse displacement field in the middle plane / line.Other modes are also present, linked to 3D effects and thus with no counterpart in the beam / plate theory, withmode shapes dominated by thickness deformations. For any mode of this second family, the displacement field in the middle plane / line has only longitudinal components, and thus no transverse component. Examples of NB modes willbe shown throughout the paper, especially in Tab. 2.Let us decompose the displacement vector X by denoting as q r , r = 1 , ..., N B the bending coordinates, and p s , s = N B +1 , ..., N the membrane coordinates: X = [ q , ..., q N B , p N B +1 , ..., p N ] T . Then, Eq. (7) can be rewritten for eachcoordinates [13,8] and involves quadratic and cubic coupling terms between the q r and the p s . We restrict ourselvesto the case of a transverse low frequency excitation, for which the external forces remain normal to the middle planeof the plate. As a consequence, the dynamics is dominated by the bending modes, which are the only ones thatreceive external excitation. In this case, Eq. (7) can be simplified. First, all quadratic α rij coefficients involving twobending coordinates i, j vanishes, in order to fulfil the symmetry of the restoring force [8,39]. In addition, one canassume that the bending coordinates, which are directly excited by the external forcing, are considered of the order A. Vizzaccaro, A. Givois, P. Longobardi, Y. Shen, J.-F. De¨u, L. Salles, C. Touz´e, O. Thomas magnitude of a small parameter ε : q r = O ( ε ), for all r ∈ { , . . . , N B } . On the other hand, NB coordinates, since theyare not directly excited and shall vibrate at a lower order of magnitude, are assumed to scale as ε : p s = O ( ε ) for all s ∈ { N B + 1 , . . . , N } . Plugging these two scaling in Eq. (7) and keeping only the leading order, one arrives to, for thebending coordinates, ∀ r = 1 , ..., N B :¨ q r + 2 ζ r ω r ˙ q r + ω r q r + N B (cid:88) i =1 N (cid:88) l = N B +1 α ril q i p l + N B (cid:88) i =1 N B (cid:88) j = i N B (cid:88) k = j β rijk q i q j q k + O ( ε ) = Q r , (12)and for the membrane coordinate, ∀ s = N B + 1 , ..., N :¨ p s + 2 ζ s ω s ˙ p s + ω s p s + N B (cid:88) i =1 N B (cid:88) j = i α sij q i q j + O ( ε ) = 0 . (13)In the above equation, one can notice that because of the transverse excitation, the second member of Eq. (13) is zero.Eqs. (12,13) approximate the dynamics of the structure in the present case of a 3D FE model. If an analyticvon K´arm´an model was used, those equations would be exact, with the terms O ( ε ) and O ( ε ) identically vanishing[8]. Those equations also show that the in-plane vibrations are quadratically coupled to the bending coordinates,while in the equations of motion of the transverse modes, only two nonlinear terms have to be taken into account: aquadratic coupling involving a product between a transverse and an in-plane coordinate, and a cubic term involvingthree transverse modes. This very specific form of equations renders the case of flat structures easier to solve thangeneral shell problems that encompass all the possible nonlinear couplings as stated in Eq. (7).2.3 Static condensation and nonlinear normal modesSince the non-bending (NB) modes have natural frequencies very large as compared to those of the directly excitedbending modes, ω s (cid:29) ω r , the dynamical part of Eqs. (13) can be neglected. The nonlinearities being more simple inthis case, one can directly express the non-bending coordinate as function of the bending ones as: p s = − N B (cid:88) i =1 N B (cid:88) j = i α sij ω s q i q j . (14)Substituting Eq. (14) into Eq. (12), one can rewrite the dynamics of the structure as a closed system involving onlybending coordinates as ¨ q r + 2 ζ r ω r ˙ q r + ω r q r + N B (cid:88) i =1 N B (cid:88) j = i N B (cid:88) k = j Γ rijk q i q j q k = Q r , (15) where the cubic Γ rijk coefficients appear. Their general expression is derived in Appendix A. If one is interested in deriving a reduced-order model for a single bending mode (say the master mode with label p ), taking into account all the other non-bending mode, then Eq. (15) can be used and the leading nonlinear cubicterm simply reduces to: Γ pppp = β pppp − N (cid:88) s = N B +1 C psppp (16)where the correction factors have been introduced and read: C psppp = α pps α spp ω s = 2( α spp ) ω s . (17) on-intrusive ROM for the dynamics of geometrically nonlinear flat structures using 3D FEM 7 These expressions show that the cubic term β pppp of the standard modal expansion must be corrected by the summationof all NB modes quadratically coupled to the master one.The quadratic coefficients α pij have some symmetry relationships, provided that the nonlinear stiffness derives froma potential [23]. In particular, the following relationship holds: α pps = 2 α spp , which leads to the second equation (17).Recalling Eq. (10), the evaluation of α spp requires only the computation of the nonlinear force f nl when the displacementalong the p -th linear mode is prescribed, whereas the calculation of the α pps coefficients for each s -th mode would requireas many evaluation of f nl as the number of non-bending modes. Nevertheless, it requires the projection of nonlinearforce onto each membrane eigenvector φ s , and thus their computation, which is costly in practice since a large numberof them is required to reach convergence (see section 3.2).From a physical perspective, following Eq. (10), coefficient α spp can be seen as the projection onto the s -th membranemode of the quadratic stiffness forces arising in the structure when a displacement along the linear p -th mode is imposed.Interestingly, this quadratic coefficient is related to a monomial q p on the s -th oscillator equation for p s . These termsare recognized as invariant-breaking terms (see e.g. [46,42]), in the sense that as soon as energy is given to the mastermode p , all s modes having these important invariant-breaking terms will no longer be vanishing. These invariant-breaking terms are responsible for the loss of invariance of the linear eigenspaces, and they are found back naturallyas correction factors when applying static condensation. They are also key in the formulation of invariant manifoldsin order to define NNMs in phase space [46,33].In parallel to the static condensation emphasised here, one can use the reduction formulae given by the normalform approach, restricting the motion to a single Nonlinear Normal Mode (NNM) [46,42,47]. In this case, the reducedorder model is directly constructed from Eqs. (7). The main advantage as compared to the above described staticcondensation is that there is no need to assume the particular structure of the equations obtained for flat structures(Eqs. (12,13)), thus generalizing the results to arches and shells. Considering only the NNM label p , the reduced-ordermodel reads:¨ q p + ω p q p + (cid:32) N (cid:88) s = N B +1 − α pps α spp ω s (cid:18) ω s − ω p ω s − ω p (cid:19) + β pppp (cid:33) q p + (cid:32) N (cid:88) s = N B +1 α pps α spp ω s (cid:18) ω s − ω p (cid:19)(cid:33) q p ˙ q p = 0 . (18)Once again, one can observe that the correction brought to the cubic term β pppp is solely given by the quadraticinvariant-breaking terms. If one has been able to compute all the quadratic α pij coefficients appearing in Eq. (18),then the model can be used to simulate the dynamics. Also, it is worth mentioning that since the NB modes havehigh frequencies, we can assume that ω s (cid:29) ω p (which is equivalent to neglect the membrane inertia). Then the termin factor of q p in Eq. (18) exactly reduces to the one obtained with static condensation in Eq. (15). On the otherhand, the term in factor of q p ˙ q p has no counterpart in static condensation, but is an order of magnitude smaller sinceit scales as 1 /ω s , so that both models are almost equivalent when a slow/fast decomposition can be assumed. Thisextends the results of [4], in which the term q p is chosen as the leading term for experimental identification purposes.A complete comparison of static condensation and nonlinear normal modes is also provided in [34] in the context ofclarifying the implicit condensation and expansion method. Finally, one can note that the formula used in Eq. (18)have been obtained thanks to a normal form approach on the conservative system [46], but they can be extended inorder to take into account the damping of the slave modes in the master coordinate ROM, hence accounting for a finer prediction of the losses [43], a feature that once again is not possible with static condensation. L = 1 m, h = 1 mm, b = 50 mm. The Young’s modulus is chosen as E = 210 GPa. This A. Vizzaccaro, A. Givois, P. Longobardi, Y. Shen, J.-F. De¨u, L. Salles, C. Touz´e, O. Thomas x yz (a) Thin beam mesh x yz (b) Thick beam mesh
Fig. 1: Mesh used in the FE computations for the first two test cases on thin and thick beams. β riii β riii α sii α ril i = r = 1 i = r = 2 i =2 , s = N B +2 i = r =2 , l = N B +2Analytic coefficient (Ac) 1.334e+03 2.128e+04 -110.0 -660.24STEP, shell elements, ν = 0 1.334e+03 2.128e+04 -110.4 -660.96Relative error with Ac (%) 0.02 % 0.005 % 0.36 % 0.11 %STEP, shell elements, ν = 0 . ν = 0 2.668e+03 4.257e+04 -109.9 -660.26Relative error with Ac (%) 100 % 100% 0.09 % 0.003 %STEP, 3D elements, ν = 0 . Table 1: Nonlinear dimensionless coefficients α ril , α sij and β rijk of the clamped beam, with non-zero and zero Poisson’sratios ( ν = 0 and ν = 0 . i = 1 , l = N B + 2). The maximum of displacement amplitudes has been fixed at h/
20 for thesecomputations.particular geometry has been chosen to be thin (the thickness to length ratio is 10 − ) to compare the results of theSTEP computation to analytic values, obtained from a beam model with Euler-Bernoulli kinematics and von K´arm´an assumptions, see [8] where these comparisons have been more fully addressed. Table 1 presents the computations ofnonlinear modal coupling coefficients obtained by the classical STEP with 3D and shell elements, compared to theanalytical values. The two meshes used here consist of four node DKQ shell elements and twenty-node brick elements(HEX20), respectively. The computations are realized with the open software Code Aster [6]. 100 elements in length,4 elements in the width have been used for both meshes, with 2 elements in the thickness for the 3D mesh. Thissufficiently refined mesh ensures that there is no convergence issue for the computations of the nonlinear coefficients.The results clearly highlights the fact that using blindly 3D elements in a STEP computation with the modal basisleads to individual values of coupling coefficients that are far from their reference, analytical values. In particular,the cubic coefficients are largely overestimated and a strong dependence to the Poisson’s ratio is found with the 3Delements: the β riii are exactly twice the expected result with 3D elements and a zero Poisson’s ratio, but they become on-intrusive ROM for the dynamics of geometrically nonlinear flat structures using 3D FEM 9 almost three times overestimated with ν = 0 .
3. On the other hand, using shell elements allows recovering the exactanalytical result if selecting ν = 0, whereas a 10 % error is found for the same STEP computation with ν = 0 .
3. Theseresults clearly demonstrate that the modal basis as input for the STEP can be used safely with 2D elements but itsextension to 3D elements is very problematic and should lead to large errors. As already noticed in the introduction, theproblem comes from the fact the one uses eigenmodeshape functions as projection basis, but not from the calculationprocedure itself.3.2 Condensation of the cubic coefficient and frequency-response curves -35 -30 -25 -20 -15 -10 -5 (a) Beam discretised with 1287 dofs -35 -30 -25 -20 -15 -10 -5 (b) Beam discretised with 5733 dofs Fig. 2: Nondimensional correction factor C s /β (Eq. (17)), associated to the first bending mode ( p = 1) and to allthe other modes ( s = 2 , . . . N ) of the thick beam testcase, over the nondimensional mode number s/N . The Poissonratio is ν = 0 . the hardening/softening behaviour of bending modes.In order to shed light on the couplings arising between the modes, another test case is chosen. It is a thick beam,with the same length L = 1 m but with a square cross section with h = b = 30 mm, as shown in Fig. 1(b). Thecross section was chosen square to be able to easily observe the 3D deformations of the cross section. The materialproperties are E = 210 GPa for the Young’s modulus and ρ = 7800 kg/m for the density. A coarse mesh of 15 HEX20elements along the axis and 2 × Γ .Fig. 2(a) shows the behaviour of the correction factor C s (defined in Eq. (17)), normalized by the cubic coefficient β , as a function of the mode number s , for the thick beam having 1287 dofs. This plot shows that the number of modes that are coupled to the first bending mode by invariant-breaking terms is very large, and uniformly distributedalong the frequency spectrum. In order to facilitate the readings, the modes for which the correction factor is below10 − have been sorted as negligible. In this family of modes that are not important, one find backs all the oddaxial modes, for symmetry reason. On the other hand, all even axial modes are strongly coupled to the first one. Themost surprising result is that if one does considers only axial modes, then only a few portion of the couplings will berevealed and taken into account. Indeed, Fig. 2(a) shows that there is a very large number of modes having very largefrequencies, and still being strongly coupled to the master bending mode.In order to check the independence of this behaviour from the mesh refinement, a second mesh of 20 elements onthe axis and 3 x 3 elements on the section has been defined on the same beam geometry. Fig. 2(b) reports a verysimilar behaviour for this second test case, where the distribution of coupled modes is uniform along the whole set ofmodes. -35 -30 -25 -20 -15 -10 -5 (a) Beam discretised with 1287 dofs -35 -30 -25 -20 -15 -10 -5 (b) Beam discretised with 5733 dofs Fig. 3: Nondimensional correction factor C s /β (Eq. (17)), associated to the first bending mode ( p = 1) and to allthe other modes ( s = 2 , . . . N ) of the thick beam testcase, over the sorted nondimensional mode number s/N , wherethe imposed order is by decreasing correction factor. The Poisson ratio is ν = 0 . − as a threshold for the significance of each contribution, a small percentage of modes, around 20%, is actually relevant. Consequently, the number of relevant modes depends on the mesh: themore refined it is, the more relevent modes are needed to reach convergence. Moreover, as seen on Fig. 2, these modesare spread over the entire spectrum, which would need the computation of all the eigenmodes of the structure, anoperation impossible in practice for a complex structure with a larger number of dofs.In order to gain insight into these coupled modes, Table 2 shows the associated mode shapes, sorted according tothe importance of their contribution in the correction factor, thus following Fig. 3(a). The table shows the first nineeigenmode shapes, recalling in the first column their number of appearance when the modes are sorted according tothe eigenfrequencies. One can observe that only one of these modes is a pure axial mode: the fourth one appearingin table 2, also being the 34 th by order of increasing frequencies. All the other ones involve important deformation inthe thickness of the beam. They are thus called thickness modes, their presence being the direct consequence of 3D on-intrusive ROM for the dynamics of geometrically nonlinear flat structures using 3D FEM 11Coupled modes ω l / (2 π )[ Hz ] C l /Γ Shape Section330 1.137e5 3.22e-1328 1.125e5 1.34e-11143 3.698e5 1.05e-134 1.043e4 9.26e-21147 3.699e5 7.01e-2370 1.202e5 4.60e-2167 7.644e4 3.35e-2172 7.653e4 2.93e-2324 1.112e5 2.89e-2 Colormap
Table 2: Order of appearance in the basis ( y -direction for the thick clamped-clamped beam. The colors scale themodulus of the displacement field and arrows display the axial displacement for axial mode 34. The Poisson ratio is ν = 0 . effects. The second column of Table 2 displays the eigenfrequencies, showing that they all are high-frequency modes.The last column shows the deformation of the section, showing the importance of thickness deformation. (a) Forced response functions at beam centre node (b) Convergence of partially corrected cubic coefficients Fig. 4: Convergence to the full model solution with the increase of coupled modes taken into account in the reducedmodel. (a) Frequency response curve of the beam at center, in the vicinity of the first bending mode eigenfrequency;case of the thick beam with 1287 dofs. Red: full model, and convergence curves with increasing number of modesretained in the truncation: 1 mode (yellow), 4 modes (green), 9 modes (blue), all modes statically condensed (dashedlight blue). Solution with one NNM in black. (b) Convergence of the evaluated corrected cubic stiffness coefficient Γ , defined in Eq. (16), with increasing number of linear modes kept in the truncation.The frequency-response curves of the thick beam in the vicinity of its first eigenfrequency is investigated in orderto illustrate how the static condensation and the NNM approach are able to retrieve the correct nonlinear behaviour.Fig. 4(a) shows the comparison of the solutions obtained by continuation, for different reduced-order models and thefull model solution. The latter has been obtained by solving all the degrees of freedom of the system with a parallelimplementation of harmonic balance method and pseudo-arc length continuation [2]; the computation of the full forcedresponse with 3 harmonics lasted approximately 36 hours. The convergence of the solution using static condensationwith an increasing number of modes to compute the correction is also shown. Despite only few modes have a very highcorrection factor C s , i.e. play a major role in the decrease of β (see Eq. (16)), it is the sum of the contributionsfrom all the coupled modes that makes the reduced model converge to the solution of the full one. In Fig. 4(b), thestrong stiffening effect coming from not having included enough coupled modes, is slowly reduced by their inclusion inthe basis; however, only the response obtained by static condensation of all coupled modes (cyan dashed) approximatesthe solution correctly (almost overlapped with the full model solution in red). On the other hand, the NNM solutionwith all the modes taken into account show also a direct convergence to the frequency-response curve. Fig. 4(b) shows the convergence of the corrected cubic coefficient Γ defined in Eq. (16) with the number ofmodes retained, i.e. the first mode plus the number of coupled modes condensed. When only one coupled mode istaken into account, the cubic coefficient β overestimates largely Γ (5.5 times): it is first explained by the classicalbending-membrane coupling effect, and secondly to Poisson effect relating to the results given in Table 1. With 9coupled modes the error on Γ is still significant (more than 60%). This strong overestimation of the cubic stiffnessvalue results in the unrealistic stiffening effect observed in the forced response. The number of coupled modes thatmust be taken into account to ensure an acceptable accuracy makes the use of STEP in its first classic formulation(i.e. with the eigenmodeshape functions as projection basis and without condensation) quite impractical: 44 modesgive a 1% error and 68 an error of 0.1%. The condensation of these modes onto the excited one becomes then a viableoption to drastically reduce the computational burden without affecting the accuracy of the solution. on-intrusive ROM for the dynamics of geometrically nonlinear flat structures using 3D FEM 13 In the previous section, we have shown that in the case of 3D elements, a strong coupling with thickness modesoccurred, thus rendering the convergence of the modal ROM particularly stringent. When one is able to compute allthe linear modes and associated coefficients, then static condensation and normal form approach can be used to finallyproduce accurate ROMs. However in most of the cases, the computation of all the linear modes, including the thicknessmodes appearing at very high frequencies, is out of reach. In this section, we investigate two alternative methods, forwhich there is no need to compute all the linear modes: static modal derivative, and a modified STEP.4.1 Static modal derivativesSections 2 and 3 were devoted to the derivation of a reduced order model from a modal point of view. In fact, a modalprojection of the quadratic nonlinear forces onto each mode φ s is required to obtain the coefficients α spp . Here we wantto introduce the concept of static modal derivatives (see [12]) because its application provides the same results as thestatic condensation of all non-bending modes, but without requiring the computation of their associated eigenvectors.The definition of modal derivatives have arisen from the recognition of the fact that in the nonlinear range, modeshapes and frequencies depend on amplitude [12]. Introducing this dependency in the eigenproblem defining the modes,one arrives at a quantity defined as the modal derivative [12,51]. Following the definition of static modal derivatives θ pr (SMD) from [13], it reads: θ pr = − K − (cid:18) ∂∂q p ∂ f nl ∂ x ( φ p q p ) (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) q p =0 · φ r . (19)When p = r a more convenient expression (in the point of view of its direct computation from a FE code) for themodal derivative θ pp writes: θ pp = − K − (cid:18) f nl ( λ φ p ) + f nl ( − λ φ p ) λ (cid:19) . (20)The equivalent general expression for θ pr , with p (cid:54) = r is provided in Appendix B. Eq. (20) shows how the SMD can beeasily computed from a set of applied static displacement, in a manner having analogies with the STEP. The term inparenthesis in Eq. (20) can be seen as the numerical second order derivatives with respect to λ of the nonlinear forcealong mode p , evaluated at the equilibrium position. This can be easily evaluated in any FE software by imposing, ona nonlinear structure, a displacement proportional to the linear p -th mode with first positive and then negative signin order to isolate the quadratic part of the nonlinear forces.The SMD can thus be seen as an added displacement vector that enriches the basis constituted by the linearmode p , in a way that takes into account the nonlinear deformation of the structure. The use of SMDs as addedvectors in the reduced order model basis is extensively documented in [40,51,35,30,13], and a complete comparison ofquadratic manifolds derived from modal derivatives with normal for theory is given in [48]. Here the focus is on therelationship between modal derivatives and non-bending modes and on the equivalence between static condensationof all non-bending modes and static condensation of modal derivatives. Given a system with geometric nonlinearities up to cubic order, the static modal derivatives related to the p -thand r -th linear modes can then be expressed in terms of the vector of quadratic coefficients α pr as: θ pr = − V Ω − α pr , (21)and the one relative to the p -th linear mode as: θ pp = − V Ω − α pp , (22)where the full matrix of eigenvectors V has been introduced. The detailed derivation of these two equations is givenin Appendix B together with the classical orthonormality properties of the matrix of eigenvectors V . By expanding Eq. (22) over all the modes and by noticing that, for a flat structure, α spp is non-zero only when s is a non-bendingmode, one obtains this important relationship (see Appendix B and [48]): θ pp = − N (cid:88) s = N B +1 φ s α spp ω s . (23)The SMD thus appears as a linear combination of coupled modes with factor − α spp /ω s ; therefore, it can be seen asa displacement field that takes into account the contribution of all non-bending modes into one equivalent vector. Toshow this property, the SMD relative to the first bending mode of the thick beam test case is depicted in Fig. 5(a). Theprojection of the SMD onto the linear modes of the system recovers Eq. (23). The contributions from the non-bendingmodes that have been identified in previous calculations, the fourth axial mode as well as various thickness modes,appear in the SMD. For each mode s , The modal amplitudes q s obtained from the projection coincides with thosefrom Eq. (23) i.e. they are equal to − α spp /ω s . These results recover and elaborate on those obtained in [13,35], whereit was shown that axial modes are contained in the SMDs of bending modes. The result is here extended to thicknessmodes and specified since the exact participation factor of each mode is made explicit. (a) Static modal derivatives (b) Non-bending modes Fig. 5: (a) Static modal derivative θ associated to the first bending mode. (b) Four first non-bending modes containedin the SMD θ , found equivalently by projecting θ on the linear mode basis, or by application of Eq. (23). Therelative modal participation factors q s of each of these modes numbered 34, 167, 328 and 330 (see also Table 1) is alsonumerically given, normalized by the total amplitude of the SMD θ ( q SMD = 1), and exactly recovers the factorsexhibited in Eq. (23).Once understood that the SMD allows gathering in a single vector the participation of all coupled modes, we wantto show how to retrieve directly, from the calculation of the SMD, the correct nonlinear behaviour of the structure,when the motion is restricted to a single master mode. In the specific case of a flat structure, Eq. (14) shows that the amplitudes of the NB modes can be explicitely related to the squared amplitude of the single master mode labeled p ,thanks to the static condensation, as: p s = − α spp ω s q p . (24)The physical displacement x ( q p ) that corresponds to the solution gathering together the master bending mode p andall its coupled NB modes can be written as: x ( q p ) = q p φ p + N (cid:88) s = N B +1 p s φ s . (25) on-intrusive ROM for the dynamics of geometrically nonlinear flat structures using 3D FEM 15 In this last equation, replacing p s by its value obtained from static condensation, Eq. (24), and then using Eq. (23)defining θ pp as a summation on the NB modes, one arrives easily at the fact that this physical displacement can beexpressed as a function of the modal coordinate plus the participation of the SMD as: x ( q p ) = q p φ p + 12 q p θ pp . (26)If one wants now to derive a reduced-order model composed of a single master coordinate (say q p here) and thatcontains the correct nonlinear behaviour, then the equation of motion would simply read:¨ q p + 2 ζ p ω p ˙ q p + ω p q p + ˜ Γ pppp q p = Q p , (27)with ˜ Γ pppp a corrected cubic coefficient. Thanks to Eq. (11b), one knows that a cubic coefficient can be found from thiscomputation, provided the imposed displacement is selected correctly. If the imposed displacement is along a linearmode, as in Eq. (11b), then one will retrieve the modal nonlinear coupling coefficient, but other choice of imposeddisplacement can be made. In particular, if one selects the one given by Eq. (26), then the cubic coefficient ˜ Γ pppp will contain the contribution of the master mode plus that of the SMD. Since Eq. (23) shows that in our particularcase (flat structure, one master mode), the SMD is completely equivalent to the static condensation of all coupledNB modes, then one will easily understand that ˜ Γ pppp = Γ pppp , the corrected cubic coefficient given in Eq. (16). Thecomplete proof of this result is provided in Appendix C.The main result, from the SMD perspective, is that if one restricts to a single master mode p , then the SMD canbe easily computed thanks to Eq. (20). Then the corrected cubic coefficient can be directly computed from: Γ pppp = φ T p ( f nl ( x ( λ )) − f nl ( x ( − λ ))) / λ , (28)where the imposed displacement is selected as in Eq. (26). This procedure allows then to find exactly the same correctedcubic coefficient of the master mode, but without resorting to the computation of all eigenvectors, as needed in thestatic condensation. It is thus a much more computationally efficient to use this methodology. Numerical examples areprovided in Section 6.1.4.2 A modified STEP for 3D elements Transverse directionMiddle line
Middle surfaceTransverse direction
Fig. 6: Examples of prescribed displacement field of the M-STEP, in the transverse direction and in the middle surface/ line of a plate / beamAs observed in the previous sections, the modal ROM associated to 3D FE discretization shows a slow convergencebecause of the couplings with very high frequency modes involving thickness deformations. Since these thickness modes are a peculiarity of the 3D model, they have no counterpart in plate or beam theory, which concentrate the kinematicaldescription on the middle plane / line. Also, using the STEP with plate / beam elements show a faster convergencesince one has to recover only the well-known coupling between bending and in-plane motions. In order to circumventthese difficulties, we propose here to modify the STEP by prescribing the displacements only on the middle line /plane of the structure, and to let free the other degrees of freedom (Fig. 6). The idea is to include automatically theeffects of NB modes, by a kind of implicit condensation of their motion, embedded into the prescribed displacement onthe middle line / plane. The obtained method is called the M-STEP, for Modified-STEP. Note that a comparable ideahas also been introduced in [15,49], but for 2D flat structures only, where only transverse motions were prescribed,leaving the other degrees of freedom free and thus building directly a condensed model.
We show in this section that it is possible to compute directly the cubic coefficients Γ rijk of Eq. (15) with a modifiedSTEP, without having to compute all the coefficients α kij and β kijl beforehand. Note that in the classical STEP, thethree components of the displacement field are prescribed to all the nodes of the FE mesh: the whole vector of unknown x is imposed and the FE code computation is just an evaluation of the internal force vector f e = f ( x ).Here, we choose to apply the STEP by prescribing the displacement field only to selected nodes and for selectedcomponents. Precisely, we denote by S the middle plane / line of the structure and n the bending direction. Tocompute Γ kppp for a given p ∈ { , . . . N B } , we choose to perform a FE computation by prescribing (i) only the transversecomponent of the bending mode φ p (ii) only on the nodes of the middle plane / line of the structure (Fig. 6). We thenprescribe the following time independent displacement to the structure: x | S , n = λ φ p | S , n , λ ∈ R . (29)where x | S , n corresponds to displacement x restricted to (i) the nodes of the FE mesh belonging to S and to (ii) itscomponents along the direction defined by n . In all the other nodes and directions, a zero forcing is prescribed. We thenuse the finite element code to solve this problem and obtain x as well as f ( x ) = f e everywhere. Since some componentsof x are not prescribed, a Newton-Raphson procedure is necessary to solve this nonlinear algebraic problem.To precise the method, we call master dofs the ones for which the displacement is prescribed, and slave dofs theother ones. The full displacement vector x and the internal forcing f can thus be decomposed as: x = (cid:20) x M x S (cid:21) , f = (cid:20) f M f S (cid:21) , (30)where the index M and S are associated to the master and slave dofs, respectively. The M-STEP consists in prescribing: (cid:40) x M = λ φ M p , f S = . (31)Then, solving the static problem f ( x ) = f e with the FE code leads to compute the internal force vector f M ( x ) on the master nodes and the displacement x S on the slave nodes. The solution of the problem then reads: x = (cid:20) λ φ M p x S (cid:21) , f = f e = (cid:20) f M (cid:21) . (32)Translated in the modal space, the above computations are close to the following situation. Prescribing via x onlythe transverse motion in the form of φ p , and because φ p is orthogonal to the other bending modes φ k , k (cid:54) = p , themodal coordinate are q p (cid:39) λ and q k (cid:39)
0. Considering the orthogonality relations associated to the stiffness matrix,this leads to assume that, for all s ∈ { , . . . N B } , s (cid:54) = p : φ T p Kx (cid:39) φ T p K ( λ φ p ) = λω p m p , φ T s Kx (cid:39) φ T s K ( λ φ p ) = 0 . (33) on-intrusive ROM for the dynamics of geometrically nonlinear flat structures using 3D FEM 17 In other words, it is assumed that the nonzero slave part x S of x is not involved in the orthogonality relations of thebending modes. Moreover, since the part of the displacement associated to longitudinal and thickness displacementsis not prescribed by x , the associated NB modal coordinates are not zero and their value depend on the nonlinearcoupling and the geometric nonlinearities. Finally, because any NB eigenmode φ s has zero displacements on the middleplane / line in the transverse direction, the modal forcing of the NB modes is exactly zero: φ s = (cid:20) φ S s (cid:21) ⇒ Q s = φ T s f e m s = 0 . (34)To summarize, one has: q p (cid:39) λ,q k (cid:39) ∀ k = 1 , . . . N, k (cid:54) = pQ s = 0 ∀ s = N B , . . . , N (35a)(35b)(35c)Error estimates of those assumptions will be introduced in section 4.2.2.Using the assumptions (35) in Eqs. (12,13) leads to: λω p + N (cid:88) s = N B +1 λα pps p s + λ β pppp = φ T p f ( x ) /m p , N (cid:88) s = N B +1 λα kps p s + λ β kppp = φ T k f ( x ) /m k , ∀ k = 1 , . . . , N B , k (cid:54) = pω s p s + λ α spp = 0 , ∀ s = N B + 1 , . . . , N (36)(37)(38)The above Eq. (38) leads to: p s = − α spp λ ω s , (39)that can be condensed into (36),(37) to give: β pppp − N (cid:88) l = N B +1 α pps α spp ω s λ = φ T p f ( x ) /m p − λω p , (40) β kppp − N (cid:88) l = N B +1 α kps α spp ω s λ = φ T k f ( x ) /m k , ∀ k (cid:54) = p. (41)One then recognizes the terms in parenthesis as the sought corrected cubic coefficients Γ ippp , defined by Eqs. (16) and(51), so that: Γ pppp = φ T p f ( x ) m p λ − ω p λ , Γ kppp = φ T k f ( x ) m k λ . (42)Consequently, all Γ ippp , i = 1 , . . . N B relies on a single nonlinear static finite elements computation, defined by Eq. (45).Other coefficients Γ kijl can be obtained in the same manner, by mixing different x M on several bending modes, followingthe classical STEP.This above described M-STEP method is a way of automatically embed in the computation the effect of all theNB modes nonlinearly coupled to the bending modes associated to Γ kijl . The essence of the method is to select theprescribed displacement x so that (i) it leaves free the degrees of freedom associated to the NB modes, so that theforcing Q s of the longitudinal modal coordinates in Eq. (38) is exactly zero and (ii) it is as orthogonal as possible tothe other bending modes than the p -th. In order to be able to quantify a priori the quality of the computation, a main idea is to check the validity of theassumptions used in the two first equations (35), which are true at first order but might deteriorate in case of anincorrect selection of master dofs. Equivalently, one can verify the orthogonality of the displacement vector x to thebending modes written in Eqs. (33). To that purpose, let us define the following errors: e pp = φ T p Kx λω p m p − , e pk = φ T k Kx λω p m p , ∀ k (cid:54) = p, (43)that should be small as compared to 1 because of Eqs. (33). If one wants to compute those errors with a FE in a nonintrusive way, Kx = f can be computed as the reaction force vector f of a linear static computation where x isprescribed to all the nodes of the FE mesh.Another check can also be performed by prescribing Eq. (31) into a linear static computation: Kx l = f e , (44)which gives: x l = (cid:20) λ φ M p x lS (cid:21) . (45)Since there are no geometrical nonlinearities, imposing φ p on the middle line / surface in the transverse directionshould result in a vector almost collinear to φ p , that is x l (cid:39) λ φ p . In particular, the slave part x lS of x l should be veryclose to φ S p , the slave part of φ p . We then define the following error: e ( ˆ y , y ) = || ˆ y − y |||| y || , (46)where || · || is the norm of vector · , and we check that e p = e ( x l , λ φ p ) is very small as compared to 1. β riii can be established: in Fig. 7, the values obtained by the direct application ofthe STEP are fitted to an heuristic law related to the Poisson ratio. In particular, the growths of the cubic coefficientsin the case of shell and 3D elements match perfectly the ratios: ρ = 11 − ν , ρ = 2(1 + ν )(1 − ν ) , (47)related to 2D and 3D constitutive laws. Indeed, the 3D constitutive law for an isotropic elastic material writes: π = E (1 + ν )(1 − ν ) (cid:2) ν tr( ε ) I + (1 − ν ) ε (cid:3) , (48)where π denotes the second Piola-Kirchhoff stress tensor, ε the Green-Lagrange strain tensor, I the identity operatorin 3D and ( E, ν ) the Young’s modulus and the Poisson ratio of the material. on-intrusive ROM for the dynamics of geometrically nonlinear flat structures using 3D FEM 190 0.1 0.2 0.3 0.4Poisson’s ratio ν β iiii / β i , A N iii , i ∈ [ , ] β β β β − ν − ν − ν
3D elementsshell elements
Fig. 7: Evolution of the ratio of the cubic coefficients β iiii compared with the analytic values β i,ANiii , with regard tothe Poisson ratio, for the first four ( i = 1 , . . .
4) bending modes of the thin beam (Fig. 1(a)) meshed with shell or3D elements, as specified on the plot. The heuristic dependences on the Poisson’s ratio are plotted with black dashedlines. All symbols are merged.Moreover, in the case of usual plate theories, a plane stress state is assumed, for which the transverse component π zz = 0 of the stress tensor is zero ( z being the direction normal to the middle plane of the plate). In this case, thein-plane counterpart of (48) reads: ˜ π = E − ν (cid:2) ν tr(˜ ε ) I + (1 − ν )˜ ε (cid:3) , (49)where ˜ π and ˜ ε denote the plane parts of π and ε , respectively, and I the identity operator in 2D.Equations (48) and (49) makes directly appear the ratios ρ and ρ , which suggests also a direct relationshipbetween the constitutive law and the results presented in Fig. 7. As a side note, the reference values, used to comparenondimensional values of cubic coefficients in Fig. 7, differ from those in Figs. 4(b) and 8(b). Indeed, in Fig. 7, thenormalising coefficient has been selected as β ,AN , the analytical value obtained from the beam theories (see e.g. [8]).As shown in the previous sections, this coefficient has to be compared with the corrected coefficient used once theconvergnce is obtained from 3D models. On the other hand, in Fig. 4(b) and 8(b), the reported values are normalizedwith respect to β , i.e. the uncorrected cubic coefficient, the one coming from direct application of STEP and shownin Table 1. This explains the different values observed between the two figures (e.g. β /β ,AN ≈ . ν = 0 . Fig. 7, whereas β /Γ = 5 . β iiii with respect to its corresponding analytical value in the case ν = 0. This leads us to investigate the couplings inthis particular case. -35 -30 -25 -20 -15 -10 -5 (a) (b) Fig. 8: Convergence without Poisson effect ( ν = 0) for the thick beam. (a) Nondimensional correction factor C s /β of all the modes over the nondimensional mode number s/N . (b) Convergence of the evaluated corrected cubic stiffnesscoefficient Γ , defined in Eq. (16), with increasing number of linear modes kept in the truncation. (c) Mode shapesand order of appearance in the basis ( ν = 0, of Figs. 3(a), 4(b) and Tab. 2. Comparing those figures shows that the number of coupled modes is much smaller in the case ν = 0 than for ν = 0 .
3: the relevant modes correspond to 5% of themodal basis if ν = 0, whereas it was 20% for ν = 0 . − in Figs 3(a)and 8(a)). Moreover, the deformations of the cross section in the case ν = 0 are purely in the bending transverse- y direction in the case ν = 0 ( y , as defined in Figs. 1 and 14, is the deformation direction of the first bending modeconsidered in all computations of the present article), whereas they were more complex (in 2D) in the case ν = 0 . ν = 0 and ν = 0 .
3. In particular, Fig. 9 shows that the deformations of the SMD cross section occurwithout distorsion: initially a square, it is deformed in a rectangle. In the case of no Poisson’s ratio, the deformationis purely in the bending y -direction, whereas the Poisson’s effects adds a slight deformation in the lateral z -direction. on-intrusive ROM for the dynamics of geometrically nonlinear flat structures using 3D FEM 21 undeformed ν=0.3ν=0 Cross section:Static modal derivatives
Fig. 9: Static modal derivative and wiew of the cross-section in the undeformed and deformed configurations, withboth ν = 0 and ν = 0 . – the deformation of the neutral axis / plane of the structure, described by classical beam / plate theories; – the 3D rotation of the cross section around the z -axis, that produces an axial deformation with a linear dependencein the thickness coordinate y ; – the 3D Poisson effect, that distorts the cross section in its two ( y and z ) directions.Then, by computing the Green-Lagrange strain tensor with this particular linear deformation, it is shown that thegeometrical nonlinearities adds two contributions to the classical von K´arm´an beam / plate model: –
3D effects that are independent of the Poisson effect , that explains a stretching in the transverse y direction,without any deformation in the lateral z direction. Those effects are a direct consequence of the 3D rotation of thecross section created by the bending; –
3D Poisson effect, that involve stretching in both the transverse y and lateral z directions.Those two geometrical effects are purely 3D and are additional to the classical membrane / bending coupling.Having in mind those observations, we can now explain those nonlinear thickness coupling. We have first to remarkthat in both cases of a beam(1D)/plate(2D) von K´arm´an model and the present 3D model, the modal expansion ofEq. (7) is exact, provided N is the number of degrees of freedom of the model. Moreover, Tab. 3 and Fig. 4(a) showthat all models converge to the same solution, which proves that the relevant modes of the basis combine themselves indifferent ways to give, at the end, the same solution. In fact, as a consequence of the above observations, the thicknessmodes are here to geometrically compensate (i) the nonlinear deformations in the transverse bending direction due tothe 3D rotation of the cross sections, shown in blue on Fig. 9 and (ii) the additional deformations of the cross sectiondue to the Poisson effect. Looking again at the deformed cross section shown in red in Fig. 9, one can understand thatat the end, the complex 3D distorsions of the cross section due to the Poisson effect (shown in Fig. 14(c)) must be fully compensated by the NB modes, which then need to be numerous. This explains the bad convergence of the modalexpansion, even worse in the case of a non zero Poisson ratio (5% of the modal basis if ν = 0, and 20% for ν = 0 . In this section, numerical examples on the three different strategies proposed in order to overcome the bias observedwhen using 3D elements, are given. In each case, the dominant cubic coefficient of a single mode is compared, usingeither the M-STEP, the static condensation of all the coupled linear modes, or the static modal derivative. Two testcases are used for the comparisons: the clamped-clamped beam of Fig. 1(b), and a clamped circular plate.
Fig. 10: First-mode displacements prescribed on the middle line of the clamped beam and the middle surface of thecircular plate. For visualisation purpose, the mesh of the plate is less fine than the one use for the computations of theROM coefficients.6.1 Application to a clamped-clamped beamComputations of the condensed coefficients Γ rijk are first performed with the three methods. For the M-STEP, the prescribed displacement is depicted in Fig. 10. The comparison made in Tab. 3 attests that the condensation with allthe eigenmodes and the first static modal derivative are quasi equivalent, whereas the M-STEP gives very close values.The relative errors are very small ( < . h/L = 0 .
03 of the beam, which is not so small to fullyverify Euler-Bernoulli assumptions.
Corrected Coefficients Γ Γ Γ Γ Γ M-StEP 2.9790e+08 -2.3555e+08 2.7964e+09 -1.9044e+09 1.9235e+10Static condensation 2.9792e+08 -2.3572e+08 2.8003e+09 -1.9083e+09 1.9288e+10Static Modal Derivative 2.9792e+08 -2.3573e+08 2.8004e+09 -1.9084e+09 1.9288e+10Analytic coefficients 2.9150e+08 -2.3151e+08 2.7082e+09 -1.8310e+09 1.8687e+10
Table 3: Corrected cubic coefficients with the three condensation methods and comparison with anaytical values. Dueto the symmetry of the mode shapes, the nonlinear coefficients which couple the first and the third modes have allnonzero values, and are therefore chosen for these computations.
In the case of the M-STEP with the displacement field prescribed on the neutral fiber, Fig. 11(a) gives the sensitivityto the prescribed displacement amplitudes of the corrected cubic coefficients Γ iiii for different modes i . It is shown thata range of validity for the displacements amplitude centered around max( λφ p ) (cid:39) h/ on-intrusive ROM for the dynamics of geometrically nonlinear flat structures using 3D FEM 23 range ofvalidity (a) Beam cross section with master nodes (b)
Fig. 11: (a) Dependence of the cubic coefficient Γ iiii with regard to the prescribed displacement amplitude, for differentmodes i , with i = 1 , , ,
4. Black dashed line : reference analytical value ( Γ i,F Eiii = Γ i,ANiii ), red dashed lines: limits ofthe range of validity. (b) Dependance of the cubic coefficient Γ with regard to the prescribed displacement amplitudefor different lines where the displacements are prescribed. Black dashed line : reference value ( Γ ,F E = Γ ,AN ). Theinset shows the location, in the cross section of the beam, of the lines in which the master displacement is prescribed.one of the lateral lines are far from the reference values. This feature will be analyzed in the following consideringerror estimators.The error estimate e pi , introduced in Section 4.2.1, is now analyzed. The criterion e pi is first computed with themaster transverse dofs prescribed in the neutral fiber, a single master mode p = 1 and different transverse modes k ∈ { , ..., } . The values presented in Fig. 12(a) show that the orthogonality is well verified, until the upper limit ofvalidity range observed on Fig. 11(a), after which the values e and e deviate from 0. For p >
1, the coefficients e pp evolve in a similar way as e , as depicted on Fig. 12(b). Tab. 4 give the numerical values of e i in the plateau ofFig. 12(a), proving that the error is the order of 10 − . Consequently, the orthogonality of the prescribed displacementfield x to the transverse modes φ i is well verified, thus validating assumptions (35a,b).Then, the same error estimate e pi is computed in the cases of a master displacement prescribed in the other lines ofthe inset of Fig. 11(b). The obtained values are given in Tab. 4. In this cases, we quantitatively confirm the observationlinked to Fig. 11(b): the orthogonality of the displacement are not verified when the master dofs are placed on a lineof the lateral surfaces of the beam. In particular, the values of the criterion e pi presented with p = 1 and i = 1 , , , in Table 4 highlight a loss of orthogonality between the first and the odd modes i = 1 , ,
5, in the case of master dofson a lateral line: indeed, the values of e , e and e deviate from 0.The second error estimate e , also introduced in Section 4.2.1 and linked to an estimate of the collinearity of theprescribed displacement x l to the master transverse modes φ , in the case of a linear computation (see Eq. (44)).This estimate confirms the above results, in particular that the displacements must preferentially be prescribed on theneutral line.A physical explanation of those effects can be deduced from the 3D displacement field of the modes. Because ofthe Poisson effect and the rotation of the sections , the displacement field on the nodes at other locations from theneutral line is not purely transverse for a bending eigenvector φ p and not zero for a NB eigenvector φ s . This explains Fig. 12: Dependence of the criterion e pi with regard to the prescribed displacement amplitude when displacements areprescribed on the neutral fiber. The values are presented on (a) for p = 1 and different modes i = 1 , , , ,
5. On (b),the values e ii are presented, also for i = 1 , , , , Neutral line Upper line Lateral line Up/Lat line e -3.496e-04 1.234e-04 -0.0210 -0.0213 e -6.993e-04 -2.183e-04 -8.136e-05 2.006e-04 e e e e Table 4: Values of the relative error e pi for different p , i , and master dofs. Nonlinear computations are performed withmax( λ φ M p )) = h/ φ s is notzero: φ s s (cid:54) = . on-intrusive ROM for the dynamics of geometrically nonlinear flat structures using 3D FEM 25 R = 0 . h = 0 . ρ = 7800 kg/m , Young’s modulus E = 210 GPa and Poisson ratio v = 0 .
3. As for the beam cases, a coarsemesh is chosen so as to compute all the modes and apply the different proposed methodologies. Consequently, 540HEX20 elements on the face and 2 HEX20 elements in the thickness were used, with a total of 1931 nodes and 4928degrees of freedom.The convergence study and appearance of thickness modes are investigated for the fundamental axisymmetricbending mode of the clamped plate, as well as the first asymmetric (1,0) mode, having one nodal line and no nodalcircle. The case of the first axisymmetric mode is awaited to share the same complexity as the beam case for symmetryreasons, but the asymmetric mode might be more difficult to achieve convergence. -35 -30 -25 -20 -15 -10 -5 (a) -35 -30 -25 -20 -15 -10 -5 (b) Fig. 13: (a) Normalised modal correction factor for the clamped circular plate, as a function of the normalised modenumber (normalization by the number of dofs), for the first asymmetric (1,0) mode of the plate. (b) The correctionfactors are now sorted by decreasing values, and two cases are shown : the case of the first asymmetic mode, corre-sponding to sorting (a), and the case of the first axisymmetric mode, showing a faster convergence. Grey points arenegligible modes in terms of coupling, magenta points are the important in-plane coupled modes while blue points arethe important thickness modes.
Fig. 13(a) shows the behaviour of the normalised modal correction factor 2( α npp /ω n ) /β pppp used in the previoussections, where p refers to the master mode (either p = 1 for the first axisymmetric mode, or p = 2 for the firstasymmetric) and n ∈ { , N } with N the number of dofs. In Fig. 13(a) only the case of the first asymmetric modeis shown for the sake of brevity (thus p = 2), but for p = 1 the trend was very similar. As for the beam, a strongcoupling with very high frequency modes is also observed. Investigating more precisely which modes are involved inthe couplings, it is found that the ones having the most important correction factor are once again thickness modes.Table 5 shows the deformed shape of the first nine modes, sorted according to their correction factor, which are thusthe most important in the coupling with the bending (1,0) mode. Two purely in-plane modes are found, in position5 and 8, and all others are thickness modes. The deformed shapes can be compared to that obtained for the beam and shown in Table 2. Indeed, the first thickness mode having the most important correction factor shows a similargeometry for both structures. Strong similarities are also observed between the second mode of the beam and mode(c) in Table 5, and the ninth mode in each case.Fig. 13(b) shows the normalized correction factor now sorted by order of decreasing values, and for the two casesof the axisymmetric fundamental mode and first asymmetric mode. It shows in particular that the convergence on thecorrect cubic coefficient is more rapidly achieved for the axisymmetric mode, where less than 10% of the modes areneeded. On the other hand, the convergence is more difficult for the first asymmetric mode. Concerning the couplingwith high-frequency modes and thickness modes for these two first bending modes, it is interesting to note that thesubset of coupled modes is almost exactly the same in the two cases, showing in particular that the coupling with thethickness modes is not very dependent on the selected bending mode. Indeed, more than 90% of the coupled modesare the same for the two cases investigated. (a1) (b) (c) (d) (e)(a2) (f) (g) (h) (i) Table 5: Mode shapes of the 9 most relevant modes coupled with the first flexural asymmetric (1,0) mode. Only twoof them are in-plane modes: (e) and (h), while all the others are thickness modes. (a2) is a side view of the top view(a1) of the first thickness mode, in order to show the strong dependence on thickness deformation.
Table 6 gathers the numerical results for the corrected cubic coefficient Γ pijk defined in Eq. (51), with p = 2 forthe first asymmetric (1,0) mode and p = 4 for the (2,0) asymmetric mode (the first bending modes being sortedby increasing frequencies, p = 1 is the fundamental axisymmetric, p = 2 , p = 4 , Γ shows a different result for the M-STEP, however the value is very smallas compared to the other ones so that this coefficient can be compared as negligible. on-intrusive ROM for the dynamics of geometrically nonlinear flat structures using 3D FEM 27Corrected cubic Coefficients Γ Γ Γ Γ Γ M-StEP 6.4763e+10 -2.8539e+09 1.5374e+11 -57.3021 3.8776e+11Static condensation 6.4762e+10 -2.8535e+09 1.5372e+11 5.3298e+03 3.8775e+11SMDs condensed 6.4762e+10 -2.8536e+09 1.5372e+11 5.3080e+03 3.8775e+11
Table 6: Corrected cubic coefficient Γ pijk , for two flexural modes, i.e. i, j, k, p ∈ [2 , In this article, a nonlinear coupling of bending modes with thickness modes of very high frequency has been exhibited,due to geometrical nonlinearities in thin flat structures. This effect adds itself to the classical longitudinal / bendingcoupling and is the cause of a very slow convergence of a reduced order model (ROM) blindly built on a modalexpansion of the nonlinear problem. It has been shown that if all eigenmodes are computed, it is possible to embedthe effect of the non-bending modes into a master bending one, thus obtaining a reduced order model composed ofonly one nonlinear Duffing oscillator. This procedure can be done either by static condensation or by a normal formreduction, equivalent to the reduction on a single nonlinear mode. Finally, two alternative methods have been proposedto overcome the problem: the use of a static modal derivative or the direct computation of the cubic coefficients byan original method, the M-STEP, inspired by the standard STEP. Those methods have been successfully verified ondedicated examples, showing equivalent results.Most of the results presented in this paper are restricted to the case of flat structures. Indeed, the specific shapeof the equations of motion (see Eqs. (12),(13)) has been used to obtain exact equivalences between different methods.One can await that the obtained results should extend to shallow curved structures. However, for more generic shellswith all the nonlinear couplings, most of the equivalences found here won’t probably hold anymore.We focused on the case of a 3D model discretized by finite elements. In section 3.1, we have shown on an examplethat some convergence problems might also appear for thin structures meshed with plate or shell elements, andhaving at least one long edge free. Our experience on thin ribbon have shown that the same kind of phenomenonappears when blindly using the STEP with the modal basis, and are again due to the loss of invariance of the modaleigenspaces. Indeed, high-frequency modes involving lateral deformations of the two free edges appeared. We alsomade computation on a circular plate with a free edge, and found circumferential modes appearing. Consequently,our finding is not restricted to 3D elements, and is completely linked to the use of the modal basis. Note that STEPcalculations can also be realized with other input functions, and this has be done in this paper e.g. in Eq. (28). The complete investigation of the analogy between this contribution, focused on 3D elements, and the problems related toplate and shell finite elements, is postponed to a future work.
Acknowledgements
The author A. Vizzaccaro is thankful to Rolls-Royce plc for the financial support. The author A. Givois isgrateful to the French Ministry of Research and Arts et Metiers Institute of Technology for their financial support through the Ph.D.Grant of the author. The author Y. Shen wishes to thank China Scholarship Council (No.201806230253) for the funding of a three-yeardoctoral position at IMSIA, ENSTA Paris and EDF Lab. The author L. Salles is thankful to Rolls-Royce plc and the EPSRC for thesupport under the Prosperity Partnership Grant ”Cornerstone: Mechanical Engineering Science to Enable Aero Propulsion Futures”,Grant Ref: EP/R004951/1. The authors thanks J. Blaˇsos for the computation of the frequency response of the full finite element modelused as reference.8 A. Vizzaccaro, A. Givois, P. Longobardi, Y. Shen, J.-F. De¨u, L. Salles, C. Touz´e, O. Thomas
Conflict of interest
The authors declare that they have no conflict of interest.
References
1. Bazant, Z., Cedolin, L.: Stability of structures. World Scientific, Singapore (2010). Third edition2. Blahoˇs, J., Vizzaccaro, A., El Haddad, F., Salles, L.: Parallel harmonic balance method for analysis of nonlinear dynamical systems.In: proc. of Turbo Expo, ASME 2020, vol. GT2020-15392, accepted (2020)3. Cadot, O., Ducceschi, M., Humbert, T., Miquel, B., Mordant, N., Josserand, C., Touz´e, C.: Wave turbulence in vibrating plates.Chapman and Hall/CRC (2016). In C. Skiadas (editor) : Handbook of Applications of Chaos theory4. Denis, V., Jossic, M., Giraud-Audine, C., Chomette, B., Renault, A., Thomas, O.: Identification of nonlinear modes using phase-locked-loop experimental continuation and normal form. Mechanical Systems and Signal Processing , 430–452 (2018). DOI10.1016/j.ymssp.2018.01.0145. Ducceschi, M., Cadot, O., Touz´e, C., Bilbao, S.: Dynamics of the wave turbulence spectrum in vibrating plates: A numericalinvestigation using a conservative finite difference scheme. Physica D , 73–85 (2014)6. Electricit´e de France: Finite element code aster rd edn. J. Wiley & Sons (2015)8. Givois, A., Grolet, A., Thomas, O., De¨u, J.F.: On the frequency response computation of geometrically nonlinear flat structuresusing reduced-order finite element models. Nonlinear Dynamics (2), 1747–1781 (2019)9. Givois, A., Tan, J.J., Touz´e, C., Thomas, O.: Backbone curves of coupled cubic oscillators in one-to-one internal resonance: bifur-cation scenario, measurements and parameter identification. Meccanica (2020). DOI 10.1007/s11012-020-01132-210. Guillot, L., Lazarus, A., Thomas, O., Vergez, C., Cochelin, B.: A purely frequency based Floquet-Hill formulation for the efficientstability computation of periodic solutions of ordinary differential systems. Journal of Computational Physics (2020). Submitted11. Holzapfel, A.G.: Nonlinear Solid Mechanics: A Continuum Approach for Engineering Science. J. Wiley & Sons (2000)12. Idelsohn, S.R., Cardona, A.: A reduction method for nonlinear structural dynamic analysis. Computer Methods in Applied Mechanicsand Engineering (3), 253 – 279 (1985)13. Jain, S., Tiso, P., Rutzmoser, J.B., Rixen, D.J.: A quadratic manifold for model order reduction of nonlinear structural dynamics.Computers and Structures , 80–94 (2017)14. Kerschen, G., Peeters, M., Golinval, J., Vakakis, A.: Non-linear normal modes, part I: a useful framework for the structuraldynamicist. Mechanical Systems and Signal Processing (1), 170–194 (2009)15. Kim, K., Khanna, V., Wang, X., , Mignolet, M.: Nonlinear reduced order modeling of flat cantilevered structures. In: Proceedingsof the 50th Structures, Structural Dynamics, and Materials Conference, AIAA Paper AIAA-2009-2492. May 4–7, Palm Springs,California (2009)16. Kim, K., Radu, A.G., Wang, X., Mignolet, M.P.: Nonlinear reduced order modeling of isotropic and functionally graded plates.International Journal of Non-Linear Mechanics , 100 – 110 (2013)17. de Klerk, D., Rixen, D.J., Voormeeren, S.: General framework for dynamic substructuring: history, review and classification oftechniques. AIAA journal (5), 1169–1181 (2008)18. Landau, L., Lifschitz, E.: Theory of Elasticity. Elsevier Butterworth Heinemann (1986). Third edition19. Lazarus, A., Thomas, O., De¨u, J.F.: Finite element reduced order models for nonlinear vibrations of piezoelectric layered beamswith applications to NEMS. Finite Elements in Analysis and Design , 35–51 (2012)20. Lewandowski, R.: Computational formulation for periodic vibration of geometrically nonlinear structures, part I: theoretical back-ground. International Journal of Solids and Structures , 1925–1947 (1997)21. Mignolet, M., Soize, C.: Stochastic reduced-order models for uncertain geometrically nonlinear dynamical systems. ComputerMethods in Applied Mechanics and Engineering , 3951–3963 (2008)22. Mignolet, M.P., Przekop, A., Rizzi, S.A., Spottswood, S.M.: A review of indirect/non-intrusive reduced order modeling of nonlineargeometric structures. Journal of Sound and Vibration , 2437–2460 (2013)23. Muravyov, A.A., Rizzi, S.A.: Determination of nonlinear stiffness with application to random vibration of geometrically nonlinearstructures. Computers and Structures (15), 1513–1523 (2003)24. Nayfeh, A.H.: Nonlinear interactions: analytical, computational and experimental methods. Wiley series in nonlinear science,New-York (2000)25. Nayfeh, A.H., Mook, D.T.: Nonlinear oscillations. John Wiley & sons, New-York (1979)26. Nayfeh, A.H., Pai, P.F.: Linear and nonlinear structural mechanics. Wiley, New-York (2004)27. Perez, R., Wang, X., Mignolet, M.P.: Prediction of displacement and stress fields of a notched panel with geometric nonlinearity byreduced order modeling. Journal of Sound and Vibration (24), 6572 – 6589 (2014)28. Perez, R., Wang, X.Q., Mignolet, M.P.: Nonintrusive Structural Dynamic Reduced Order Modeling for Large Deformations: En-hancements for Complex Structures. Journal of Computational and Nonlinear Dynamics (3) (2014)on-intrusive ROM for the dynamics of geometrically nonlinear flat structures using 3D FEM 2929. Rizzi, S.A., Przekop, A.: System identification-guided basis selection for reduced-order nonlinear response analysis. Journal of Soundand Vibration (3), 467 – 485 (2008)30. Rutzmoser, J.B., Rixen, D.J., Tiso, P., Jain, S.: Generalization of quadratic manifolds for reduced order modeling of nonlinearstructural dynamics. Computers and Structures , 196–209 (2017)31. Salen¸con, J.: Handbook of Continuum Mechanics. Springer (2001)32. S´en´echal, A.: R´eduction de vibrations de structure complexe par shunts pi´ezo´electriques. application aux turbomachines. Ph.D.thesis, Conservatoire National des Arts et M´etiers, Paris, France (2011). URL
33. Shaw, S.W., Pierre, C.: Non-linear normal modes and invariant manifolds. Journal of Sound and Vibration (1), 170–173 (1991)34. Shen, Y., B´ereux, N., Frangi, A., Touz´e, C.: Reduced order models for geometrically nonlinear structures: comparison of implicitcondensation and nonlinear normal modes. Computers and Structures p. submitted (2020)35. Sombroek, C.S.M., Tiso, P., Renson, L., Kerschen, G.: Numerical computation of nonlinear normal modes in a modal derivativesubspace. Computers & Structures , 34 – 46 (2018)36. Szemplinska-Stupnicka, W.: The behavior of nonlinear vibrating systems, vol. I: fundamental concepts and methods. Applicationto single-degree-of-freedom systems. Kluwer Academic Publishers, Dordrecht (1990)37. Thomas, O., Bilbao, S.: Geometrically nonlinear flexural vibrations of plates: In-plane boundary conditions and some symmetryproperties. Journal of Sound and Vibration (3), 569–590 (2008)38. Thomas, O., S´en´echal, A., De¨u, J.F.: Hardening/softening behaviour and reduced order modelling of nonlinear vibrations of rotatingcantilever beams. Nonlinear dynamics (2), 1293–1318 (2016). DOI 10.1007/s11071-016-2965-039. Thomas, O., Touz´e, C., Chaigne, A.: Non-linear vibrations of free-edge thin spherical shells: modal interaction rules and 1:1:2internal resonance. International Journal of Solids and Structures (11-12), 3339–3373 (2005)40. Tiso, P., Jansen, E., Abdalla, M.: Reduction method for finite element nonlinear dynamic analyses of shells. AIAA Journal (10),2295 – 2304 (2011)41. Tiso, P., Rixen, D.J.: Discrete empirical interpolation method for finite element structural dynamics. In: Topics in NonlinearDynamics, Volume 1, pp. 203–212. Springer New York, New York, NY (2013)42. Touz´e, C.: Normal form theory and nonlinear normal modes: theoretical settings and applications. In: G. Kerschen (ed.) ModalAnalysis of nonlinear Mechanical Systems, pp. 75–160. Springer Series CISM courses and lectures, vol. 555, New York, NY (2014)43. Touz´e, C., Amabili, M.: Non-linear normal modes for damped geometrically non-linear systems: application to reduced-order mod-eling of harmonically forced structures. Journal of Sound and Vibration (4-5), 958–981 (2006)44. Touz´e, C., Bilbao, S., Cadot, O.: Transition scenario to turbulence in thin vibrating plates. Journal of Sound and Vibration (2),412–433 (2012)45. Touz´e, C., Thomas, O.: Non-linear behaviour of free-edge shallow spherical shells: effect of the geometry. International Journal ofNon-linear Mechanics (5), 678–692 (2006)46. Touz´e, C., Thomas, O., Chaigne, A.: Hardening/softening behaviour in non-linear oscillations of structural systems using non-linearnormal modes. Journal of Sound and Vibration (1-2), 77–101 (2004)47. Touz´e, C., Vidrascu, M., Chapelle, D.: Direct finite element computation of non-linear modal coupling coefficients for reduced-ordershell models. Computational Mechanics (2), 567–580 (2014)48. Vizzaccaro, A., Salles, L., Touz´e, C.: Comparison of nonlinear mappings for reduced-order modeling of vibrating structures: normalform theory and quadratic manifold method with modal derivatives. Nonlinear Dynamics p. submitted (2020)49. Wang, X., O’Hara, P., Mignolet, M., Hollkamp, J.: Reduced-order modeling with local enrichment for the nonlinear geometricresponse of a cracked panel. AIAA Journal (1), 421–436 (2019)50. Wang, X., Phlipot, G.P., Perez, R.A., Mignolet, M.P.: Locally enhanced reduced order modeling for the nonlinear geometric responseof structures with defects. International Journal of Non-Linear Mechanics , 1 – 7 (2018)51. Weeger, O., Wever, U., Simeon, B.: On the use of modal derivatives for nonlinear model order reduction. International Journal forNumerical Methods in Engineering (13), 1579–1602 (2016)52. Yuan, J., El-Haddad, F., Salles, L., Wong, C.: Numerical assessment of reduced order modeling techniques for dynamic analysis ofjointed structures with contact nonlinearities. Journal of Engineering for Gas Turbines and Power (3) (2019) A Definition of cubic nonlinear terms
A particular feature of the manipulations of nonlinear coupling coefficients lies in the fact that they are related to monoms havingsymmetry relationships. For example, a cubic coefficient β pijj is related to the monom q i q j , which is also the case of β pjij and β pjji .Consequently many formulations use upper triangular forms for quadratic and cubic coefficients α pij (assuming j ≥ i ) and β pijk (with k ≥ j ≥ i ). However, direct calculations produce coefficients that have not this ordering property built-in. The formula given here allowsone to get from to another formulation.The corrected cubic coefficients in Eq. (15) can be derived by replacing the value of p s from Eq. (14) into the quadratic term ofEq. (12): N B (cid:88) i =1 N (cid:88) s = N B +1 α ris q i p s = N B (cid:88) i =1 N B (cid:88) j =1 N B (cid:88) k = j N (cid:88) s = N B +1 − α ris α sjk ω s q i q j q k i , j , and k consistent with the ones of the cubic term of Eq. (12), i.e. havingthe sum over j from i to N B , reads: N B (cid:88) i =1 N B (cid:88) j =1 N B (cid:88) k = j N (cid:88) s = N B +1 − α ris α sjk ω s q i q j q k = N B (cid:88) i =1 N B (cid:88) j = i +1 N B (cid:88) k = j +1 N (cid:88) s = N B +1 (cid:18) − α ris α sjk ω s − α rjs α sik ω s − α rks α sij ω s (cid:19) q i q j q k + N B (cid:88) i =1 N B (cid:88) k = i +1 N (cid:88) s = N B +1 (cid:18) − α ris α sik ω s − α rks α sii ω s (cid:19) q i q k + N B (cid:88) i =1 N B (cid:88) k = i +1 N (cid:88) s = N B +1 (cid:18) − α ris α skk ω s − α rks α sik ω s (cid:19) q i q k + N B (cid:88) i =1 N (cid:88) s = N B +1 (cid:18) − α ris α sii ω s (cid:19) q i . (50)This term can be rewritten in a more compact form by defining the correction factor: C rsijk = α ris α sjk ω s + α rjs α sik ω s + α rks α sij ω s i < j < kα ris α sik ω s + α rks α sii ω s i = j < kα ris α skk ω s + α rks α sik ω s i < j = kα ris α sii ω s i = j = k leading to: N B (cid:88) i =1 N B (cid:88) j =1 N B (cid:88) k = j N (cid:88) s = N B +1 − α ris α sjk ω s q i q j q k = N B (cid:88) i =1 N B (cid:88) j = i N B (cid:88) k = j N (cid:88) s = N B +1 −C rsijk q i q j q k now with summation indexes consistent with the ones of Eq. (12).Finally, the quadratic and cubic nonlinear terms in Eq. (12) read: N B (cid:88) i =1 N (cid:88) s = N B +1 α ris q i p s + N B (cid:88) i =1 N B (cid:88) j = i N B (cid:88) k = j β rijk q i q j q k = N B (cid:88) i =1 N B (cid:88) j = i N B (cid:88) k = j β rijk − N (cid:88) s = N B +1 C rsijk q i q j q k . and the corrected cubic coefficient in Eq. (15): Γ rijk = β rijk − N (cid:88) s = N B +1 C rsijk (51) B Expression of static modal derivatives in terms of quadratic coupling coefficients
Given the general equation of a system with quadratic and cubic nonlinearities in physical coordinates: M ¨ x + C ˙ x + Kx + f nl ( x ) = f e , (52)the nonlinear force can be written in terms of nonlinear tensors as: f nl ( x ) = Axx + Bxxx , (53)on-intrusive ROM for the dynamics of geometrically nonlinear flat structures using 3D FEM 31where we used the compact tensor notation also employed in [13,30]. In order to explicit the notation, the products Axx and
Bxxx are here given with explicit indicial notation:
Axx = N (cid:88) i =1 N (cid:88) j =1 A ij x j x i , Bxxx = N (cid:88) i =1 N (cid:88) j =1 N (cid:88) k =1 B ijk x i x j x k . The most inner product defined above coincides with a matrix product performed on the last index of the tensors.In order to derive Eq. (23), we first express the i -th column of nonlinear stiffness matrix as: (cid:18) ∂ f nl ∂ x (cid:19) i = N (cid:88) j =1 A ij x j + N (cid:88) j =1 A ji x j + N (cid:88) j =1 N (cid:88) k =1 B ijk x j x k + N (cid:88) j =1 N (cid:88) k =1 B jik x j x k + N (cid:88) j =1 N (cid:88) k =1 B jki x j x k . (55)Exploiting the symmetry of the tensors A and B which implies that A ij = A ji and similarly B ijk = B jik = B jki , that stems from thefact that geometric nonlinear forces can be derived from a potential (see [23]), the nonlinear stiffness matrix can be written in compactform as: ∂ f nl ∂ x = 2 Ax + 3 Bxx . (56)The nonlinear stiffness matrix evaluated along mode p reads: ∂ f nl ∂ x ( Φ p q p ) = 2 q p AΦ p + 3 q p BΦ p Φ p , (57)and its derivatives with respect to q p evaluated at q p = 0 reads: (cid:18) ∂∂q p ∂ f nl ∂ x ( Φ p q p ) (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) q p =0 = 2 AΦ p . (58)Hence the static modal derivatives defined in Eq. (19) leads to: θ pr = − K − AΦ p Φ r , (59)and, when p = r , to: θ pp = − K − AΦ p Φ p . (60)We now want to relate these expressions to the quadratic modal coupling coefficients α spr obtained from the STEP method, as wellas expliciting directly how to compute the modal derivative from specific static evaluation of the internal force vector. When p = r ,Eq. (11b) reads: α spp = Φ T s (cid:32) f nl ( q p Φ p ) + f nl ( − q p Φ p )2 q p (cid:33) /m s . (61)When p (cid:54) = r , the STEP method needs to call two static evaluations with q p Φ p and q s Φ s so that the expression writes: α spr = Φ T s (cid:18) f nl ( q p Φ p + q r Φ r ) + f nl ( − q p Φ p − q r Φ r ) − f nl ( q p Φ p ) − f nl ( − q p Φ p ) − f nl ( q r Φ r ) − f nl ( − q r Φ r )2 q p q r (cid:19) /m s . (62)Let us now assume that the eigenvectors are mass normalised so that m s = 1 , ∀ s , using Eq. (53) we can find the relation between α spp and A as: α spp = Φ T s AΦ p Φ p . (63)Similarly for the coefficients α spr with p < r : α spr = 2 Φ T s AΦ p Φ r , (64)where the factor 2 appears for symmetry reasons and is related to the usual problem of representing the polynomial monoms bycounting them separately or not. Indeed, in the usual polynomial representation, the coefficients α spr with p < r contains both values of Φ T s AΦ p Φ r and Φ T s AΦ r Φ p whereas the coefficients α srp are set to zero. We can now observe that each coefficient α spr can be seen asthe row of a vector α pr whose compact expression would be: α pr = 2 V T AΦ p Φ r , (65)2 A. Vizzaccaro, A. Givois, P. Longobardi, Y. Shen, J.-F. De¨u, L. Salles, C. Touz´e, O. Thomaswhere the full matrix of eigenvector V has been introduced. In the same line, the expression for the vector α pp when p = r has thesame shape but without the factor 2: α pp = V T AΦ p Φ p . (66)Introducing Eqs. (65)-(66) respectively in Eqs. (59)-(60), and using the relationships of the quadratic coupling coefficients from theSTEP method, Eqs. (62)-(61), one obtains the important formulas allowing one to compute directly the modal derivative from staticFE calculations. The expression for θ pp is given in the main text as Eq. (20) while the formula for θ pr reads: θ pr = − K − (cid:18) f nl ( λ ( Φ p + Φ r )) + f nl ( − λ ( Φ p + Φ r )) − f nl ( λ Φ p ) − f nl ( − λ Φ p ) − f nl ( λ Φ r ) − f nl ( − λ Φ r )2 λ (cid:19) (67)We now want to relate more closely the modal derivative to the modal representation and derive an explicit expression showing thatthe SMD gathers the contributions of all coupled modes. For that purpose, one needs to express the usual orthonormality conditionsshared by the matrix of eigenvectors V , assumed to be mass normalized: V T MV = I , V T KV = Ω , where Ω is the diagonal matrix containing the squared eigenfrequencies.Recalling eq. (59) we can now express the static modal derivatives in terms of α pr : θ pr = − K − V − T α pr , (68)and for p = r : θ pp = − K − V − T α pp . (69)Moreover, using the orthogonality conditions: K − V − T = V Ω − , (70)we can express the static modal derivatives in the general case as: θ pr = − V Ω − α pr , (71) θ pp = − V Ω − α pp . (72)In the particular case of a flat structure, the equations of motions have a simple structure recalled in Eqs. (12)-(13), so that α spr is nonzeroonly for the non-bending modes. Thanks to this simplification, one can express the static modal derivatives as a linear combination ofthe non-bending modes only as: θ pr = − N (cid:88) s = N B +1 Φ s α spr ω s , (73) θ pp = − N (cid:88) s = N B +1 Φ s α spp ω s . (74) C Corrected cubic coefficient obtained from static modal derivatives
This appendix aims at demonstrating that coefficient ˜ Γ pppp introduced in Eq. (27), i.e. by using an imposed displacement composed ofthe master mode plus a quadratic part containing the SMD, is exactly equal to that given in Eq. (16), and obtained thanks to the staticcondensation. For that purpose, let us recall that the imposed displacement in the first case reads: x ( q p ) = q p Φ p + 12 q p θ pp . (75)The cubic coefficient ˜ Γ pppp can be computed by using the general formula from the STEP method, Eq. (11b), by replacing the imposeddisplacement by the one given in Eq. (75). Consequently, one arrives at:˜ Γ pppp = Φ T p (cid:18) f nl ( λ Φ p + 12 λ θ pp ) − f nl ( − λ Φ p + 12 λ θ pp ) (cid:19) / λ . (76)on-intrusive ROM for the dynamics of geometrically nonlinear flat structures using 3D FEM 33Using the explicit expression of the nonlinear forces from Eq. (53): f nl ( λ Φ p + 12 λ θ pp ) = λ AΦ p Φ p + λ (cid:18) AΦ p θ pp + 12 Aθ pp Φ p + BΦ p Φ p Φ p (cid:19) + O ( λ ) (77) f nl ( − λ Φ p + 12 λ θ pp ) = λ AΦ p Φ p − λ (cid:18) AΦ p θ pp + 12 Aθ pp Φ p + BΦ p Φ p Φ p (cid:19) + O ( λ ) . (78)Thus the difference between the above nonlinear forces reads: f nl ( λ Φ p + 12 λ θ pp ) − f nl ( − λ Φ p + 12 λ θ pp ) = 2 λ (cid:18) AΦ p θ pp + 12 Aθ pp Φ p + BΦ p Φ p Φ p (cid:19) + O ( λ ) (79)and the coefficient ˜ Γ pppp , neglecting high order terms, reads:˜ Γ pppp = 12 Φ T p AΦ p θ pp + 12 Φ T p Aθ pp Φ p + Φ T p BΦ p Φ p Φ p (80)Following a similar argument as the one of Eq. (61) it is possible to show that the last term on the right hand side of Eq. (80) is theuncorrected cubic coefficient: β pppp = Φ T p BΦ p Φ p Φ p (81)As regards to the first two terms on the right hand side of Eq. (80), using the link between SMD and non-bending modes of Eq. (23)valid in case of a flat structure, they can be written as:12 Φ T p AΦ p θ pp = − N (cid:88) s = N B +1 Φ T p AΦ p Φ s α spp ω s (82)12 Φ T p Aθ pp Φ p = − N (cid:88) s = N B +1 Φ T p AΦ s Φ p α spp ω s (83)where the index s spans over all the non-bending modes. Using the symmetry of the tensor A once again: Φ T p AΦ s Φ p + Φ T p AΦ p Φ s = 2 Φ T p AΦ p Φ s (84)and in light of Eq. (64): 2 Φ T p AΦ p Φ s = α pps (85)because p > s for each non-bending mode.It is now possible to express the corrected cubic coefficient as:˜ Γ pppp = β pppp − N (cid:88) s = N B +1 α pps α spp ω s = Γ pppp (86)and recover the same expression of the corrected coefficient given in Eq. (16) obtained with static condensation of all non-bendingmodes.It is worth mentioning that the demonstration could have been done with another pathway. Indeed, the displacement introduced inEq. (75) follows the general strategy proposed in [13,30], consisting of using SMD to build a quadratic manifold approach in order todefine a reduced-order models thanks to a nonlinear mapping between physical and reduced coordinates. However, the SMD can alsobe used more simply, by considering θ pp as an enrichment of the modal basis, composed here of the single master eigenvector Φ p . Onecould then derive a two-dofs reduced-order model, by projecting the general equations of motion onto these two vectors. Then, sincethe motions associated to the SMD are linked to NB modes having high frequencies, one can neglect their inertia and proceed to thestatic condensation of the part coming from the static modal derivative. By doing so, one would show again that the corrected cubiccoefficient is still equal to Γ pppp . Hence the static condensation of the SMD (a single vector) is thus strictly equivalent, in our simplifiedcase of a flat structure (and considering only one master mode), to the static condensation of all coupled NB modes (including thicknessmodes).4 A. Vizzaccaro, A. Givois, P. Longobardi, Y. Shen, J.-F. De¨u, L. Salles, C. Touz´e, O. Thomas (a) S M b σ xx e x e z e y x y (b) Magni fi ed with a factor 20 yz(c) Fig. 14: Analytical solution of a beam in pure bending. (a) Stress state σ xx in a cross section S ; (b,c) reference (ingrey) and deformed (in blue) configurations of the beam. The Poisson effect has been magnified with a factor 20 inthe cross section view (c). D Analytical solution for the pure bending of a beam
We consider the linear elastic solution of a beam of rectangular cross section under pure bending (see Fig. 14). The beam is thussubjected to a uniform bending moment M b = M e z . The material of the beam is assumed linear elastic with Young’s modulus E and Poisson’s ratio ν . The orthonormal frame ( e x , e y , e z ) is used, with e x colinear to the middle axis of the beam, e y the direction ofbending and e z = e x ∧ e y with ∧ the vector product. The local equilibrium of the beam is exactly verified by the following axial stressstate, linear through the thickness of the beam:div σ = ⇒ σ = σ xx , σ xx = − αy, (87)where div is the divergence operator, σ is the stress tensor and α is a constant. The bending moment writes M = αI where I is thesecond moment of inertia of the beam. The linear strain tensor ε verifies: ε = 12 (cid:16) ∇ U + ∇ T U (cid:17) = 1 + νE σ − νE tr σ I ⇒ (cid:15) = − α y/E − νy/E
00 0 − νy/E (88)As exposed in [31], the following displacement field verifies exactly the above equations: U = αxy e x − α (cid:2) x + ν ( y − z ) (cid:3) e y − ναyz e z . (89)This displacement is composed of three parts: – the transverse ( e y ) component, proportional to x , which is the standard transverse displacement of the neutral line due to thebending. This term is the one directly computed in a beam theory; – the axial ( e x ) component, which is the 3D linearized rotation of the cross section around vector e z , also due to bending. In a beamtheory, it is the consequence of the Euler-Bernoulli kinematics; – two additional terms in the transverse ( e y ) and lateral ( e z ) directions, proportional to ν and thus directly linked to the Poissoneffect. These terms are responsible of the distortion of the cross section, as seen in Fig 14(c).on-intrusive ROM for the dynamics of geometrically nonlinear flat structures using 3D FEM 35To see the effect of the geometrical nonlinearities, we use the displacement field (89) to compute the nonlinear Green-Lagrange straintensor. We obtain: γ = 12 (cid:16) ∇ U + ∇ T U + ∇ T U ∇ U (cid:17) (90)= ε + α x (cid:124) (cid:123)(cid:122) (cid:125) γ + α y xy xy x
00 0 0 (cid:124) (cid:123)(cid:122) (cid:125) γ + να xy − xzxy ν ( y + z ) 0 − xz ν ( y + z ) (cid:124) (cid:123)(cid:122) (cid:125) γ (91)The above equation shows that, in addition to the linear part ε , the nonlinear part has three type of components: – the first nonlinear term γ is purely axial, with the γ xx term α x / v ( x ) = α/ x denotesthe transverse displacement of the neutral fiber, the nonlinear terms added by the von K´arman theory in γ xx is v (cid:48) ( x ) / α x / – the second nonlinear term γ comes from purely 3D effects independent of the Poisson effect , that add (i) a stretch in the axial e x direction, proportional to y ; (ii) a positive and homogeneous stretch in the transverse e y direction (proportional to x ); (iii) atransverse shear; – the third nonlinear term γ gathers the effects of the 3D Poisson effect. It involves stretching in both the transverse e y and lateral e zz