Nonlinear vibrations of large structures with uncertain parameters

Abstract

The effects of uncertainties on the nonlinear dynamics of complex structures remain poorly mastered and most methods deal with the linear case. This article deals with a model of a large and complex structure with uncertain parameters for the nonlinear dynamic case, and the reduction in the model discretized by the finite element method is obtained by reducing the degrees of freedom in the numerical model. This is achieved by the development of the unknown displacement vector on the basis of the eigenmodes; a particular attention is paid to the calculation of the nonlinear stiffness coefficients of the model. The method combines the stochastic finite element methods with a modal reduction class based on sub-structuring the component mode synthesis method. The reference method is the Monte Carlo simulation which consists in making several simulations for different values of the uncertain parameters. The simulation of complex and nonlinear structures is costly in terms of memory and computation time. To solve this problem, the perturbation method combined with the component mode synthesis reduction method significantly reduces the computational cost by preserving the physical content of the original structure. The numerical integration by the Newmark schema is used; the first statistical moments (mean and variance) of the nonlinear dynamic response are computed. Numerical simulations illustrate the accuracy and effectiveness of the proposed methodology.

Introduction

In a robust design process, the determination of the variability of the nonlinear dynamic response of a large and complex structure is essential. Uncertainties come from the tolerances of manufacturing, the boundary conditions, and the external excitations. These structures are largely used in the fields of aerospace, automotive, civil engineering, and so on.

The uncertainty of the physical parameters, nonlinearity, and complexity of the structure require the development of a complete mathematical approach for predicting the dynamic behavior variability.

In addition, several methods have been developed in the literature to take account of uncertain parameters in the nonlinear dynamic response.

The reference method is the Monte Carlo simulation. ¹ This method allows a statistical evaluation based on a large number of deterministic analyses by considering different values of uncertain parameters. However, it requires the generation of big size samples and then generates a prohibitory time computing. This approach is thus very costly.

Also, perturbation methods are widely used to calculate the first moments (mean and standard deviation) of dynamic response whose uncertain variables vary slightly. These techniques are based on the Taylor series development of the response around its mean. This allows the direct determination of the variability of the response according to the physical parameters (mechanical and geometrical) randomly.

Indeed, perturbation methods based on a development in Taylor series of second order ² and Neumann expansion method ³ are generally efficient. Another development in the first order ⁴ gives a similar result to the previous developments with a reduced time computing. Recently, in the field of thermal conduction with uncertain parameters,^5–7 combination of a probabilistic study based on the perturbation method with a numerical resolution based on the finite difference method is applied. ⁵ A prediction of the temperature field with random and fuzzy parameters of the properties of materials is studied; a numerical technique called fuzzy stochastic finite element method based on the combination of the perturbation theory with the moment methods is adopted. ⁶

Furthermore, another form of development is a polynomial chaos expansion (PCE).^8,9 The stochastic solution may be expanded in terms of the polynomial chaos basis whose elements are obtained from orthogonal polynomial. ¹⁰ The properties of this polynomial basis are used to generate a system of deterministic equations. The resolution of this system is used to determine the variability of the response.

However, nonlinear dynamics, the large number of degrees of freedom (DOFs) due to the mesh of a large structure, and higher order developing for modeling uncertainty induced a considerable increase in deterministic equations.

One way to solve this problem is the reduction by component mode synthesis (CMS) method proposed in the literature.^11–18 This method allows condensing the large number of DOFs into a small number using the generalized coordinates. FA Lulf et al. ¹⁹ proposed a comparative study of different bases for reduction in nonlinear dynamic structures. Thus, in the CMS method, the overall structure is divided into sub-structures, each of which is analyzed independently in order to obtain the corresponding solution. These solutions are combined to obtain the overall solution of the structure by imposing constraints on the interfaces. The different methods are classified according to CMS interface: fixed interface, ¹¹ free interface,^12,13 or hybrid interface.^14,16

Recently, D Sarsri et al.^20,21 developed an approach coupling CMS reduction method and developing uncertainty by a PCE to calculate the frequency transfer functions and response temporal for linear stochastic structures. J Sinou et al. ²² proposed for simple structures, requiring no reduction, a technique taking into account the uncertainties in nonlinear models by combining the method of harmonic balance method (HBM) and developing uncertainty by a PCE. This method is based on a formulation of nonlinear dynamic problem in which the physical parameters, nonlinear forces, and the excitation force are considered randomly.

In another work, D Sarsri and L Azrar ²³ used the CMS method coupled with the perturbation method to calculate the stochastic modes of large finite element models with uncertain parameters for the linear problems.

This work is an extension to the nonlinear problems. The aim is to estimate the stochastic nonlinear dynamic response for a large structure with a minimum computational cost. To do this, we develop a methodological approach for calculating the temporal response of a large structure with uncertain parameters. This approach is based on coupling of the perturbation method and the reduction (CMS) method. First, we develop the nonlinear dynamic equations considering geometrical nonlinearity. The resolution of the nonlinear dynamic problem by the finite element method is adopted. Then, the temporal integration by Newmark scheme is developed. Second, we take completely the random phenomena using the perturbation method. The method of stochastic finite element is used. Various types of CMS interface method are used to optimally reduce the model size. The first moments of the nonlinear dynamic response of the reduced system are compared with the entire system. Several numerical simulations have shown the accuracy and efficiency of procedures and methodologies developed.

Modeling of nonlinearity

The discretization by the finite element method of a linear structure gives the following matrix system

$$[M]\{\stackrel{··}{u}\}+[C]\{\stackrel{·}{u}\}+[K]\{u\}=\left\{F_{\mathrm{ext}}\right\}$$

(1)

where $[M]$ is the mass matrix, $[C]$ is the damping matrix, [K] is the stiffness matrix, $\{u\}$ is the displacement vector, and $\{F_{\mathrm{ext}}\}$ is the external force vector. In this system there may be nonlinear phenomena due to the following:

• Intermittent contact or friction called nonlinearity of contact modeled by nonlinear force

  $$\left\{F_{\mathrm{nl}}(\{u\})\right\}=\{\begin{array}{c}{k}_1\left\{u\right\}\mathrm{if}\{u\}>0\hfill \\ \multicolumn{1}{c}{k_2\left\{u\right\}\mathrm{if}\{u\}<0}\end{array}$$

• Large displacement for thin structures called geometric nonlinearity modeled by nonlinear force $\{F_{\mathrm{nl}}(\{u\})\}=[K_{\mathrm{nl}}]\{u{\}}^p$ , p essentially quadratic or cubic.

The nonlinear system modeling the nonlinear dynamics of a structure discretized by the finite element method is given by

$$[M]\{\stackrel{··}{u}\}+[C]\{\stackrel{·}{u}\}+[K]\{u\}+\left\{F_{\mathrm{nl}}(\{u\})\right\}=\left\{F_{\mathrm{ext}}\right\}$$

(2)

where $\{F_{\mathrm{nl}}(\{u\})\}$ is the nonlinear force vector. In this study, we used a cubic geometric nonlinearity.

CMS reduction method

The CMS method consists in using simultaneously a sub-structuring technique and a reduction method. The large and complex structure is partitioned into sub-structures. Each sub-structure is represented by a reduced basis composed of the normal modes and the interface modes. We present the theoretical bases of the CMS method. Initially, the eigenmodes and the interface static deformations are given for each sub-structure. Then, the overall system is projected on these bases taking into account the interface couplings between the sub-structures, after the reduced system is solved. Finally, the complete system solution is reconstituted.

The finite element model of the entire structure is partitioned into N sub-structures SS(i) (i = 1, …, N). The equations of motion for each nonlinear sub-structure SS(i) are

$$[M{]}^i\{\stackrel{··}{u}{\}}^i+[C{]}^i\{\stackrel{·}{u}{\}}^i+[K{]}^i\{u{\}}^i+{\left\{F_{\mathrm{nl}}\right\}}^i={\left\{F_e\right\}}^i$$

(3)

where $[M{]}^i$ , $[C{]}^i$ , and $[K{]}^i$ are, respectively, the mass matrix, the damping matrix, and the stiffness matrix for sub-structures SS(i).

The displacement vector $\{u{\}}^i$ is partitioned into a vector $\{u_j{\}}^i$ , called interface DOF and $\{u_{\mathrm{in}}{\}}^i$ is the vector of internal DOF

$$\{u{\}}^i={\left\{\begin{array}{c}{u}_j\hfill \\ \multicolumn{1}{c}{u_{\mathrm{in}}}\end{array}\right\}}^i$$

(4)

The external force vector $\{F_e{\}}^i$ is composed of vectors $\{F_{\mathrm{ej}}{\}}^i$ and $\{F_{\mathrm{ee}}{\}}^i$ , called interface force and external applied force

$${\left\{F_e\right\}}^i={\left\{F_{\mathrm{ej}}\right\}}^i+{\left\{F_{\mathrm{ee}}\right\}}^i$$

(5)

The nonlinear force vector $\{F_{\mathrm{nl}}{\}}^i$ is composed of vectors $\{F_{\mathrm{nlj}}{\}}^i$ and $\{F_{\mathrm{nle}}{\}}^i$ , called interface force and external nonlinear force

$${\left\{F_{\mathrm{nl}}\right\}}^i={\left\{F_{\mathrm{nlj}}\right\}}^i+{\left\{F_{\mathrm{nle}}\right\}}^i$$

(6)

In the CMS methods, the physical displacements of the sub-structure SS(i) are expressed as a linear combination of the sub-structure modes. After some algebraic transformations, a set of Ritz vectors Q is obtained and the displacement vector of each sub-structure can be expressed as

$${\left\{u\right\}}^{(i)}={\left[Q\right]}^{(i)}\left\{\begin{array}{c}{u}_j^{(i)}\hfill \\ \multicolumn{1}{c}{{\eta }_p^{(i)}}\end{array}\right\}={\left[Q\right]}^{(i)}{\left\{u_c\right\}}^{(i)}$$

(7)

where ${\eta }_p^{(i)}$ are the generalized coordinates. The matrix $[Q{]}^{(i)}$ is defined according to the method of sub-structuring used (fixed or free interface ¹⁵ ).

• 1. Free interface method

In the free interface method, the displacements of each sub-structure are expressed as

$$\left\{u^{(i)}\right\}=\left[{\Phi }^{(i)}\right]\left\{{\eta }^{(i)}\right\}+\left[{\psi }_r^{(i)}\right]\left\{{\xi }_r^{(i)}\right\}+\left[{\psi }_a^{(i)}\right]\left\{{\xi }_a^{(i)}\right\}$$

(8)

where $[{\Phi }^{(i)}]$ is the matrix of truncated undamped normal modes of the sub-structure $S{S}^{(i)}$ with a free interface as boundary condition. $[{\psi }_r^{(i)}]$ is the matrix of rigid body modes for an unconstrained sub-structure $S{S}^{(i)}$ with a free interface. $[{\psi }_a^{(i)}]$ is the matrix of attachment modes associated with the interface, which are the static deformation shapes of $S{S}^{(i)}$ obtained by applying successively a unit force to one coordinate of the interface

$$\left[{\psi }_a^{(i)}\right]=\left[G^{(i)}\right]\left\{F_j^{(i)}\right\}$$

(9)

where $\{F_j^{(i)}\}={\left\{\begin{array}{c}{I}_j\hfill \\ \multicolumn{1}{c}{0}\end{array}\right\}}^{(i)}$ and $[G^{(i)}]$ is the residual flexibility matrix.

The expression of $[G^{(i)}]$ depends on the nature of the problem. If the sub-structure is statically determined (i.e. no rigid body modes), then

$$\left[G^{(i)}\right]={\left[K^{(i)}\right]}^{-1}$$

else

$$\left[G^{(i)}\right]={}^t\left[A^{(i)}\right]{\left[K_c^{(i)}\right]}^{-1}\left[A^{(i)}\right]$$

with

$$\left[A^{(i)}\right]=\left[I\right]-\left[{\phi }^{(i)}\right]{}^t\left[{\phi }^{(i)}\right]$$

and

$${}^t\left[{\phi }^{(i)}\right]\left[M^{(i)}\right]\left[{\phi }^{(i)}\right]=\left[I\right]$$

where $[I]$ is an unit matrix. $[{\phi }^{(i)}]$ is a matrix of rigid modes. $[K_c^{(i)}]$ is a stiffness matrix obtained by fixing arbitrary DOF to make the structure isostatic and replacing the corresponding part of the initial stiffness matrix by zero.

To preserve the interface DOF, the following partition is used

$$\begin{array}{c}\multicolumn{1}{c}{\left[{\Phi }^{(i)}\right]=\left[\begin{array}{c}{\Phi }_j^{(i)}\hfill \\ \multicolumn{1}{c}{{\Phi }_{\mathrm{in}}^{(i)}}\end{array}\right]}\\ \multicolumn{1}{c}{\left[{\psi }_r^{(i)}\right]=\left[\begin{array}{c}{\psi }_{rj}^{(i)}\hfill \\ \multicolumn{1}{c}{{\psi }_{r\mathrm{in}}^{(i)}}\end{array}\right]}\\ \multicolumn{1}{c}{\left[{\psi }_a^{(i)}\right]=\left[\begin{array}{c}{\psi }_{aj}^{(i)}\hfill \\ \multicolumn{1}{c}{{\psi }_{a\mathrm{in}}^{(i)}}\end{array}\right]}\end{array}$$

Using this partition, one obtains

$$\begin{array}{c}\multicolumn{1}{c}{\left\{{{\xi }_a}^{\left(i\right)}\right\}={\left[{\psi }_{a_j}^{(i)}\right]}^{-1}\left\{u_j^{(i)}\right\}-{\left[{\psi }_{a_j}^{(i)}\right]}^{-1}}\\ \left[{\Phi }_j^{(i)}\right]\left\{{\eta }^{(i)}\right\}-{\left[{\psi }_{a_j}^{(i)}\right]}^{-1}\left[{\psi }_{r_j}^{(i)(i)}\right]\left\{{\xi }_r^{(i)}\right\}\hfill \end{array}$$

(10)

The matrix Q is then given by

$${\left[Q\right]}^i=\left[\begin{array}{c}{\psi }_a^{(i)}{\psi }_{a_j}^{{(i)}^{-1}}{\psi }_r^{(i)}-{\psi }_a^{(i)}{\psi }_{a_j}^{{(i)}^{-1}}{\psi }_{r_j}^{(i)}\hfill \\ \multicolumn{1}{c}{{\Phi }^{(i)}-{\psi }_a^{(i)}{\psi }_{a_j}^{{(i)}^{-1}}{\Phi }_j^{(i)}}\end{array}\right]$$

(11)

• 2. Fixed interface method

In the fixed interface method, the displacements of each sub-structure are expressed as

$$\left\{u^{(i)}\right\}=\left[{\Phi }^{(i)}\right]\left\{{\eta }^{(i)}\right\}+\left[{\psi }_c^{(i)}\right]\left\{u_j^{(i)}\right\}$$

(12)

The matrix $[Q{]}^i$ is given by

$${\left[Q\right]}^i=\left[{\psi }_c^{(i)}{\Phi }^{(i)}\right]$$

(13)

where $[{\Phi }^{(i)}]$ is the matrix of truncated undamped normal modes of the sub-structure $S{S}^{(i)}$ with a fixed interface as boundary condition. $[{\psi }_c^{(i)}]$ is the matrix of the constrained mode associated with the interface, which is the static deformation shape of $S{S}^{(i)}$ obtained by imposing successively a unit displacement on one interface while holding the remaining interface coordinates fixed.

The conservation of interface DOF allows assembling these matrices as in the ordinary finite element methods. Let us denote the vector of independent displacements of the assembled structure by $\{u_c\}$

$$\left\{u_c\right\}=\left\{\begin{array}{c}\begin{array}{c}{\eta }_p^{(1)}\hfill \\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{{\eta }_p^{(N)}}\end{array}\hfill \\ \multicolumn{1}{c}{u_j}\end{array}\right\}$$

(14)

The compatibility of interface displacements of the assembled structure is obtained by writing for each sub-structure SS(i) the following relation

$${\left\{u_c\right\}}^i={\left[\beta \right]}^i\left\{u_c\right\}$$

(15)

where $[\beta {]}^i$ is the matrix of localization or of geometrical connectivity of the SS(i) sub-structure. It makes possible to locate the DOF of each sub-structure SS(i) in the global DOF of the assembled structure. They are the Boolean matrices whose elements are 0 or 1.

A transformation matrix can be defined for each sub-structure SS(i) by

$${\left[Z\right]}^i={\left[Q\right]}^i{\left[\beta \right]}^i$$

(16)

where $[Q{]}^i$ is given by the considered CMS method.

The displacement vector $\{u{\}}^{(i)}$ is then given by

$${\left\{u\right\}}^{\left(i\right)}={\left[Z\right]}^i\left\{u_c\right\}$$

(17)

Inserting equation (17) into equation (1) and multiply on the right by ${}^t{[Z]}^i$ , using the sum for all sub-structures, the following equation is obtained

$$\begin{array}{c}\multicolumn{1}{c}{\left[M_c\right]\left\{{\stackrel{··}{u}}_c\right\}+\left[C_c\right]\left\{{\stackrel{·}{u}}_c\right\}+\left[K_c\right]\left\{u_c\right\}}\\ +\sum _{i=1}^N{}^t{\left[Z\right]}^i({\left\{F_{\mathrm{nlj}}\right\}}^i+{\left\{F_{\mathrm{nle}}\right\}}^i)\hfill \\ =\sum _{i=1}^N{}^t{\left[Z\right]}^i({\left\{F_{\mathrm{ej}}\right\}}^i+{\left\{F_{\mathrm{ee}}\right\}}^i)\hfill \end{array}$$

(18)

where

$$\begin{array}{c}\multicolumn{1}{c}{\left[M_c\right]=\sum _{i=1}^N{}^t{\left[Z\right]}^i{\left[M\right]}^i{\left[Z\right]}^i}\\ \multicolumn{1}{c}{\left[C_c\right]=\sum _{i=1}^N{}^t{\left[Z\right]}^i{\left[C\right]}^i{\left[Z\right]}^i}\\ \multicolumn{1}{c}{\left[K_c\right]=\sum _{i=1}^N{}^t{\left[Z\right]}^i{\left[K\right]}^i{\left[Z\right]}^i}\end{array}$$

(19)

Using the interface DOF compatibility of displacements, it can easily be shown that

$$\sum _{i=1}^N{}^t{\left[Z\right]}^i{\left\{F_{\mathrm{ej}}\right\}}^i=0$$

(20)

Finally, the reduced equation of motion can be written as follows

$$\begin{array}{c}\multicolumn{1}{c}{\left[M_c\right]\left\{{\stackrel{··}{u}}_c\right\}+\left[C_c\right]\left\{{\stackrel{·}{u}}_c\right\}+\left[K_c\right]\left\{u_c\right\}}\\ +\sum _{i=1}^N{}^t{\left[Z\right]}^i({\left\{F_{\mathrm{nlj}}\right\}}^i+{\left\{F_{\mathrm{nle}}\right\}}^i)\hfill \\ =\sum _{i=1}^N{}^t{\left[Z\right]}^i{\left\{F_{\mathrm{ee}}\right\}}^i\hfill \end{array}$$

(21)

Stochastic perturbation method

The perturbation method is largely employed in the field of the stochastic finite elements. It is based on an approximation the random variables by their development in Taylor series around their average value. These developments are truncated at the second order. The perturbation method must obey the conditions of existence and validity, in particular the reduced field of variation of the random variables. We present an extension of this method for the nonlinear dynamic systems with uncertain parameters.

Let us assume that for each sub-structure, the mass matrix $[M{]}^i$ , the dumped matrix $[C{]}^i$ , the stiffness matrix $[K{]}^i$ , and the external force vector $\{F_e{\}}^i$ are related to a vector of the random variables ${\vartheta }_i$ (i = 1, …, I). Thus, the condensed mass $[M_c]$ , dumped $[C_c]$ , and stiffness $[K_c]$ matrices are related to same vector of the random variables.

In the time domain, the resulting reduced stochastic differential system equation (21) has to be solved. The first two moments of time response (average and variance) will be calculated using the second-order perturbation method.

One defines the vector of the average parameters $\overline{{\vartheta }_i}$ and the quantity $d{\vartheta }_i={\vartheta }_i-\overline{{\vartheta }_i}$ . All the matrices and vector in equation (21) are random and are expanded through second-order Taylor series as follows

$$\begin{array}{c}\left[M_c\right]={\left[M_c\right]}^0+{\left[M_c\right]}^{n}d{\vartheta }_n+{\left[M_c\right]}^{\mathrm{np}}d{\vartheta }_{n}d{\vartheta }_p\hfill \\ \left[C_c\right]={\left[C_c\right]}^0+{\left[C_c\right]}^{n}d{\vartheta }_n+{\left[C_c\right]}^{\mathrm{np}}d{\vartheta }_{n}d{\vartheta }_p\hfill \\ \left[K_c\right]={\left[K_c\right]}^0+{\left[K_c\right]}^{n}d{\vartheta }_n+{\left[K_c\right]}^{\mathrm{np}}d{\vartheta }_{n}d{\vartheta }_p\hfill \\ \left[F_e\right]={\left[F_e\right]}^0+{\left[F_e\right]}^{n}d{\vartheta }_n+{\left[F_e\right]}^{\mathrm{np}}d{\vartheta }_{n}d{\vartheta }_p\hfill \end{array}$$

(22)

where [.]⁰, [.] ⁿ , and [.] ^np are deterministic matrices corresponding to the zero-, the first-, and the second-order partial derivatives with respect to the random parameter ${\vartheta }_i$ and given by

$$\begin{array}{c}{\left[X\right]}^0={\left[X\right]\left(\vartheta \right)|}_{\vartheta =\overline{\vartheta }}\hfill \\ {\left[X\right]}^n={\frac{\partial \left[X\right]\left(\vartheta \right)}{\partial {\vartheta }_n}|}_{\vartheta =\overline{\vartheta }}\hfill \\ {\left[X\right]}^{\mathrm{np}}=\frac{1}{2}{\frac{{\partial }^2\left[X\right]\left(\vartheta \right)}{\partial {\vartheta }_n\partial {\vartheta }_p}|}_{\vartheta =\overline{\vartheta }}\hfill \end{array}$$

Indicial notations are used, with indices n, p running over the sequence 1, 2, …, I as well as the repeated indices summation.

For structures with small uncertainties, one can assume that the transformation matrix $[Z]$ is deterministic. The zero-, first-, and second-order derivatives of the condensed mass $[M_c]$ , dumped $[C_c]$ , and stiffness $[K_c]$ matrices are given by

$$\begin{array}{c}{\left[X_c\right]}^0=\sum _{i=1}^N{}^t{\left[Z\right]}^{(i)}{\left[X\right]}^{0(i)}{\left[Z\right]}^{(i)}\hfill \\ {\left[X_c\right]}^n=\sum _{i=1}^N{}^t{\left[Z\right]}^{(i)}{\left[X\right]}^{n(i)}{\left[Z\right]}^{(i)}\hfill \\ {\left[X_c\right]}^{\mathrm{np}}=\sum _{i=1}^N{}^t{\left[Z\right]}^{(i)}{\left[X\right]}^{\mathrm{np}(i)}{\left[Z\right]}^{(i)}\hfill \end{array}$$

(23)

where $[X_c]$ can take the symbol $[M_c]$ , $[C_c]$ , or $[K_c]$ .

The unknown vector displacement, velocity, and acceleration are also developed through Taylor series as follows

$$\begin{array}{c}\left\{u_c\right\}={\left\{u_c\right\}}^0+{\left\{u_c\right\}}^{n}d{\vartheta }_n+{\left\{u_c\right\}}^{\mathrm{np}}d{\vartheta }_{n}d{\vartheta }_p\hfill \\ \left\{{\stackrel{·}{u}}_c\right\}={\left\{{\stackrel{·}{u}}_c\right\}}^0+{\left\{{\stackrel{·}{u}}_c\right\}}^{n}d{\vartheta }_n+{\left\{{\stackrel{·}{u}}_c\right\}}^{\mathrm{np}}d{\vartheta }_{n}d{\vartheta }_p\hfill \\ \left\{{\stackrel{··}{u}}_c\right\}={\left\{{\stackrel{··}{u}}_c\right\}}^0+{\left\{{\stackrel{··}{u}}_c\right\}}^{n}d{\vartheta }_n+{\left\{{\stackrel{··}{u}}_c\right\}}^{\mathrm{np}}d{\vartheta }_{n}d{\vartheta }_p\hfill \end{array}$$

(24)

And the vector nonlinear force is given as

$$\begin{array}{c}\left\{F_{\mathrm{nl}}(\{u\})\right\}=\left\{F_{\mathrm{nl}}({\{u\}}^0+{\{u\}}^{n}d{\vartheta }_n+{\{u\}}^{\mathrm{np}}d{\vartheta }_{n}d{\vartheta }_p)\right\}\hfill \\ \left\{F_{\mathrm{nl}}(\{u\})\right\}=\left\{F_{\mathrm{nl}}\left({\{u\}}^0\right)\right\}\hfill \\ +\left\{\frac{\partial F_{\mathrm{nl}}}{\partial u}({\{u\}}^{n}d{\vartheta }_n+{\{u\}}^{\mathrm{np}}d{\vartheta }_{n}d{\vartheta }_p)\right\}\hfill \end{array}$$

(25)

Then

$$\begin{array}{c}{\left\{F_{\mathrm{nl}}\right\}}^0=\left\{F_{\mathrm{nl}}\left({\{u\}}^0\right)\right\}\hfill \\ {\left\{F_{\mathrm{nl}}\right\}}^n=\frac{\partial F_{\mathrm{nl}}}{\partial u}\left({\{u\}}^n\right)\hfill \\ {\left\{F_{\mathrm{nl}}\right\}}^{\mathrm{np}}=\frac{\partial F_{\mathrm{nl}}}{\partial u}\left({\{u\}}^{\mathrm{np}}\right)\hfill \end{array}$$

The condensed nonlinear forces are

$$\begin{array}{c}{\left\{F_{\mathrm{nl}}\right\}}_c^0=\sum _{i=1}^N{}^t{\left[Z\right]}^{\left(i\right)}{\left\{F_{\mathrm{nl}}\right\}}^{0(i)}\hfill \\ {\left\{F_{\mathrm{nl}}\right\}}_c^n=\sum _{i=1}^N{}^t{\left[Z\right]}^{\left(i\right)}{\left\{F_{\mathrm{nl}}\right\}}^{n(i)}\hfill \\ {\left\{F_{\mathrm{nl}}\right\}}_c^{\mathrm{np}}=\sum _{i=1}^N{}^t{\left[Z\right]}^{\left(i\right)}{\left\{F_{\mathrm{nl}}\right\}}^{\mathrm{np}(i)}\hfill \end{array}$$

(26)

The condensed external forces are

$$\begin{array}{c}{\left\{F_e\right\}}_c^0=\sum _{i=1}^N{}^t{\left[Z\right]}^{(i)}{\left\{F_e\right\}}^{0(i)}\hfill \\ {\left\{F_e\right\}}_c^n=\sum _{i=1}^N{}^t{\left[Z\right]}^{(i)}{\left\{F_e\right\}}^{n(i)}\hfill \\ {\left\{F_e\right\}}_c^{\mathrm{np}}=\sum _{i=1}^N{}^t{\left[Z\right]}^{(i)}{\left\{F_e\right\}}^{\mathrm{np}(i)}\hfill \end{array}$$

(27)

Substituting these developments into equation (21) and writing the terms of same order on gets the following differential systems:

• Zero-order equation

  $${\left[M_c\right]}^0{\left\{{\stackrel{··}{u}}_c\right\}}^0+{\left[C_c\right]}^0{\left\{{\stackrel{·}{u}}_c\right\}}^0+{\left[K_c\right]}^0{\left\{u_c\right\}}^0+{\left\{F_{\mathrm{nl}}\right\}}_c^0={\left\{F_e\right\}}_c^0$$
  (28)

• First-order equation

  $$\begin{array}{c}\multicolumn{1}{c}{{\left[M_c\right]}^0{\left\{{\stackrel{··}{u}}_c\right\}}^n+{\left[C_c\right]}^0{\left\{{\stackrel{·}{u}}_c\right\}}^n+{\left[K_c\right]}^0{\left\{u_c\right\}}^n+{\left[M_c\right]}^n{\left\{{\stackrel{··}{u}}_c\right\}}^0}\\ +{\left[C_c\right]}^n{\left\{{\stackrel{·}{u}}_c\right\}}^0+{\left[K_c\right]}^n{\left\{u_c\right\}}^0+{\left\{F_{\mathrm{nl}}\right\}}_c^n={\left\{F_e\right\}}_c^n\hfill \end{array}$$
  (29)

• Second-order equation

  $$\begin{array}{c}\multicolumn{1}{c}{{\left[M_c\right]}^0{\left\{{\stackrel{··}{u}}_c\right\}}^{\mathrm{np}}+{\left[C_c\right]}^0{\left\{{\stackrel{·}{u}}_c\right\}}^{\mathrm{np}}+{\left[K_c\right]}^0{\left\{u_c\right\}}^{\mathrm{np}}+{\left[M_c\right]}^{\mathrm{np}}{\left\{{\stackrel{··}{u}}_c\right\}}^0}\\ +{\left[C_c\right]}^{\mathrm{np}}{\left\{{\stackrel{·}{u}}_c\right\}}^0+{\left[K_c\right]}^{\mathrm{np}}{\left\{u_c\right\}}^0+2{\left[M_c\right]}^n{\left\{{\stackrel{··}{u}}_c\right\}}^p\hfill \\ +2{\left[C_c\right]}^n{\left\{{\stackrel{·}{u}}_c\right\}}^p+2{\left[K_c\right]}^n{\left\{u_c\right\}}^p\hfill \\ +{\left\{F_{\mathrm{nl}}\right\}}_c^{\mathrm{np}}={\left\{F_e\right\}}_c^{\mathrm{np}}\hfill \end{array}$$
  (30)

Stochastic temporal response

The temporal response from time 0 to time T of equations (28)–(30) is required. The time T is subdivided into n intervals $\Delta t=T/n$ , and the numerical solution is obtained at times $t_r=r.\Delta t$ $r\in \mathrm{IN}$ and $0\le r\le n$ Assume that the solutions at times t are known and that the solution at time $(t+\Delta t)$ is required next. According to the Newmark method, the following assumption is used at time $(t+\Delta t)$

$$\begin{array}{c}\multicolumn{1}{c}{{\{{\stackrel{··}{u}}_c(t+\Delta t)\}}^x=\frac{{\{u_c(t+\Delta t)\}}^x-{\{u_c(t)\}}^x}{\alpha {(\Delta t)}^2}}\\ -\frac{{\{{\stackrel{·}{u}}_c(t)\}}^x}{\alpha (\Delta t)}-\frac{(1-2\alpha ){\{{\stackrel{··}{u}}_c(t)\}}^x}{2\alpha }\hfill \\ \multicolumn{1}{c}{{\{{\stackrel{·}{u}}_c(t+\Delta t)\}}^x=\frac{\delta ({\{u_c(t+\Delta t)\}}^x-{\{u_c(t)\}}^x)}{\alpha (\Delta t)}}\\ -(1-\frac{\delta }{\alpha }){\{{\stackrel{·}{u}}_c(t)\}}^x-(\Delta t)(1-\frac{\delta }{2\alpha }){\{{\stackrel{··}{u}}_c(t)\}}^x\hfill \end{array}$$

(31)

where x can take 0, n, and np values.

In which, the two parameters $\alpha$ and $\delta$ verify $\delta \ge 1/2$ and $\alpha \ge (\delta +0.5)/4$ in order to get accurate and stable solution. The following notations are used

$$\begin{array}{c}{a}_0=\frac{\delta }{\alpha {(\Delta t)}^2}{a}_1=\frac{\delta }{\alpha (\Delta t)}\hfill \\ a_2=\frac{1}{\alpha \left(\Delta t\right)}{a}_3=\frac{1}{2\alpha }-1\hfill \\ a_4=\frac{\delta }{\alpha }-1a_5=\frac{(\Delta t)}{2}(\frac{\delta }{\alpha }-1)\hfill \\ a_6=(\Delta t)(1-\delta )a_7=(\Delta t)\delta \hfill \end{array}$$

Based on these notations, the following equations are resulted:

• Zero-order equation

$${\left[K_{\mathrm{eq}}\right]}^0{\left\{u_c(t+\Delta t)\right\}}^0={\left\{F_{\mathrm{eq}}\right\}}^0$$

(32)

with

$$\begin{array}{c}{\left[K_{\mathrm{eq}}\right]}^0={\left[K_c\right]}^0+{\left[K_{\mathrm{nlc}}\right]}^0+a_0{\left[M_c\right]}^0+a_1{\left[C_c\right]}^0\hfill \\ {\left[K_{\mathrm{nlc}}\right]}^0=\sum _{i=1}^N{}^t{\left[Z\right]}^{(i)}{\left[K_{\mathrm{nl}}\right]}^{0(i)}{\left[Z\right]}^{(i)}\hfill \\ {\left[K_{\mathrm{nl}}\right]}^{0(i)}={\frac{\partial {\left\{F_{\mathrm{nl}}\right\}}^i}{\partial u}|}_{\left\{u\right\}={\left\{u\right\}}^{0(i)}}\hfill \\ {\left\{F_{\mathrm{eq}}\right\}}^0={\{F_e(t+\Delta t)\}}_c^0+{\left[M_c\right]}^0\hfill \\ (a_0{\{u_c(t)\}}^0+a_2{\{{\stackrel{·}{u}}_c(t)\}}^0+a_3{\{{\stackrel{··}{u}}_c(t)\}}^0)\hfill \\ +{\left[C_c\right]}^0(a_1{\{u_c(t)\}}^0+a_4{\{{\stackrel{·}{u}}_c(t)\}}^0+a_5{\{{\stackrel{··}{u}}_c(t)\}}^0)\hfill \end{array}$$

• First-order equation

  $${\left[K_{\mathrm{eq}}\right]}^n{\{u_c(t+\Delta t)\}}^n={\left\{F_{\mathrm{eq}}\right\}}^n$$
  (33)

with

$$\begin{array}{c}{\left[K_{eq}\right]}^n={\left[K_c\right]}^n+{\left[K_{nlc}\right]}^n+a_0{\left[M_c\right]}^n+a_1{\left[C_c\right]}^n\hfill \\ {\left[K_{nlc}\right]}^n={\displaystyle \sum }_{i=1}^N{}^t{\left[Z\right]}^{(i)}{\left[K_{nl}\right]}^{n(i)}{\left[Z\right]}^{(i)}\hfill \\ {\left[K_{nl}\right]}^{n(i)}={\frac{\partial {\left\{F_{nl}\right\}}^i}{\partial u}|}_{\left\{u\right\}={\left\{u\right\}}^{n(i)}}\hfill \\ {\left\{F_{eq}\right\}}^n={\left\{F_e(t+\Delta t)\right\}}_c^n+{\left[M_c\right]}^0\left(a_0{\left\{u_c(t)\right\}}^n+a_2{\left\{{\stackrel{·}{u}}_c(t)\right\}}^n+a_3{\left\{{\stackrel{··}{u}}_c(t)\right\}}^n\right)\hfill \\ +{\left[C_c\right]}^0\left(a_1{\left\{u_c(t)\right\}}^n+a_4{\left\{{\stackrel{·}{u}}_c(t)\right\}}^n+a_5{\left\{{\stackrel{··}{u}}_c(t)\right\}}^n\right)-{\left[M_c\right]}^n{\left\{{\stackrel{··}{u}}_c(t+\Delta t)\right\}}^0+{\left[C_c\right]}^n{\left\{{\stackrel{·}{u}}_c(t+\Delta t)\right\}}^0\hfill \\ +{\left[K_c\right]}^n{\left\{u_c(t+\Delta t)\right\}}^0+{\left\{F_{nl}(t+\Delta t)\right\}}_c^n\hfill \end{array}$$

• Second-order equation

$${\left[K_{\mathrm{eq}}\right]}^{\mathrm{np}}{\{u_c(t+\Delta t)\}}^{\mathrm{np}}={\left\{F_{\mathrm{eq}}\right\}}^{\mathrm{np}}$$

(34)

with

$$\begin{array}{c}{\left[K_{\mathrm{eq}}\right]}^{\mathrm{np}}={\left[K_c\right]}^{\mathrm{np}}+{\left[K_{\mathrm{nlc}}\right]}^{\mathrm{np}}+a_0{\left[M_c\right]}^{\mathrm{np}}+a_1{\left[C_c\right]}^{\mathrm{np}}\hfill \\ {\left[K_{\mathrm{nlc}}\right]}^{\mathrm{np}}=\sum _{i=1}^N{}^t{\left[Z\right]}^{(i)}{\left[K_{\mathrm{nl}}\right]}^{\mathrm{np}(i)}{\left[Z\right]}^{(i)}\hfill \\ {\left[K_{\mathrm{nl}}\right]}^{\mathrm{np}(i)}={\frac{\partial {\left\{F_{\mathrm{nl}}\right\}}^i}{\partial u}|}_{\left\{u\right\}={\left\{u\right\}}^{\mathrm{np}(i)}}\hfill \\ {\left\{F_{\mathrm{eq}}\right\}}^{\mathrm{np}}={\{F_e(t+\Delta t)\}}_c^{\mathrm{np}}+{\left[M_c\right]}^0\hfill \\ (a_0{\{u_c(t)\}}^{\mathrm{np}}+a_2{\{{\stackrel{·}{u}}_c(t)\}}^{\mathrm{np}}+a_3{\{{\stackrel{··}{u}}_c(t)\}}^{\mathrm{np}})\hfill \\ +{\left[C_c\right]}^0(a_1{\{u_c(t)\}}^{\mathrm{np}}+a_4{\{{\stackrel{·}{u}}_c(t)\}}^{\mathrm{np}}+a_5{\{{\stackrel{··}{u}}_c(t)\}}^{\mathrm{np}})·\hfill \\ (a_1{\{u_c(t)\}}^{\mathrm{np}}+a_4{\{{\stackrel{·}{u}}_c(t)\}}^{\mathrm{np}}+a_5{\{{\stackrel{··}{u}}_c(t)\}}^{\mathrm{np}})\hfill \\ -{\left[M_c\right]}^{\mathrm{np}}{\{{\stackrel{··}{u}}_c(t+\Delta t)\}}^0-{\left[C_c\right]}^{\mathrm{np}}{\{{\stackrel{·}{u}}_c(t+\Delta t)\}}^0\hfill \\ -{\left[K_c\right]}^{\mathrm{np}}{\{u_c(t+\Delta t)\}}^0\hfill \\ -{\{F_{\mathrm{nl}}(t+\Delta t)\}}_c^{\mathrm{np}}-2{\left[M_c\right]}^n{\{{\stackrel{··}{u}}_c(t+\Delta t)\}}^p\hfill \\ -2{\left[C_c\right]}^n{\{{\stackrel{·}{u}}_c(t+\Delta t)\}}^p\hfill \\ -2{\left[K_c\right]}^n{\{u_c(t+\Delta t)\}}^p-{\{F_{\mathrm{nl}}(t+\Delta t)\}}_c^n\hfill \end{array}$$

The solution of the problem is obtained by successively solving the following algebraic equations

$$\begin{array}{c}{\{u_c(t+\Delta t)\}}^0={\left({\left[K_{\mathrm{eq}}\right]}^0\right)}^{-1}{\left\{F_{\mathrm{eq}}\right\}}^0\hfill \\ {\{u_c(t+\Delta t)\}}^n={\left({\left[K_{\mathrm{eq}}\right]}^n\right)}^{-1}{\left\{F_{\mathrm{eq}}\right\}}^n\hfill \\ {\{u_c(t+\Delta t)\}}^{\mathrm{np}}={\left({\left[K_{\mathrm{eq}}\right]}^{\mathrm{np}}\right)}^{-1}{\left\{F_{\mathrm{eq}}\right\}}^{\mathrm{np}}\hfill \end{array}$$

(35)

The derivative of the physical displacements of each sub-structure is then obtained by

$$\begin{array}{c}{\{u(t+\Delta t)\}}^{0(i)}={\left[Z\right]}^{(i)}{\{u_c(t+\Delta t)\}}^0\hfill \\ {\{u(t+\Delta t)\}}^{n(i)}={\left[Z\right]}^{(i)}{\{u_c(t+\Delta t)\}}^n\hfill \\ {\{u(t+\Delta t)\}}^{\mathrm{np}(i)}={\left[Z\right]}^{(i)}{\{u_c(t+\Delta t)\}}^{\mathrm{np}}\hfill \end{array}$$

(36)

The mean and the variance values of displacement $u(t+\Delta t)$ are given by

$$\begin{array}{c}E\left(\{u(t+\Delta t)\}\right)={\{u(t+\Delta t)\}}^0\hfill \\ +\frac{1}{2}{\{u(t+\Delta t)\}}^{\mathrm{np}}\mathrm{cov}({\vartheta }_n,{\vartheta }_p)\hfill \\ \mathrm{Cov}\left(\{u(t+\Delta t)\}\right)\hfill \\ =\left({\{u(t+\Delta t)\}}^n{\{u(t+\Delta t)\}}^p\right)\mathrm{cov}({\vartheta }_n,{\vartheta }_p)\hfill \\ +\left({\{u(t+\Delta t)\}}^n{\{u(t+\Delta t)\}}^{\mathrm{pq}}\right)E({\vartheta }_n,{\vartheta }_p,{\vartheta }_q)\hfill \\ +\left({\{u(t+\Delta t)\}}^{\mathrm{np}}{\{u(t+\Delta t)\}}^{\mathrm{ql}}\right)E({\vartheta }_n,{\vartheta }_p)E({\vartheta }_q,{\vartheta }_l)\hfill \\ +E({\vartheta }_n,{\vartheta }_q)E({\vartheta }_p,{\vartheta }_l)\hfill \end{array}$$

(37)

If the correlation between the random parameters is not considered, equation (37) can be simplified as

$$\begin{array}{c}E\left(\{u(t+\Delta t)\}\right)={\{u(t+\Delta t)\}}^0\hfill \\ +\frac{1}{2}{\{u(t+\Delta t)\}}^{\mathrm{nn}}\mathrm{var}\left({\vartheta }_n\right)\hfill \\ \mathrm{Var}\left(\{u(t+\Delta t)\}\right)={\left({\{u(t+\Delta t)\}}^n\right)}^2\mathrm{Var}\left({\vartheta }_n\right)\hfill \\ +\frac{1}{4}{\left({\{u(t+\Delta t)\}}^{\mathrm{np}}\right)}^2\mathrm{Var}\left({\vartheta }_n\right)\mathrm{Var}\left({\vartheta }_p\right)\hfill \end{array}$$

(38)

Numerical example

For nonlinear discrete systems with stochastic parameters, some benchmark tests are elaborated to demonstrate the efficiency of the methodological approach. The presented method can be applied to continuous or discrete systems. In this article, we restrict ourselves to show the applicability and effectiveness of these methods for the dynamic analysis of nonlinear discrete systems with N DOFs. A nonlinear dynamic system consisting of 20 masses connected by 21 springs nonlinearly is shown in Figure 1. This structure will be divided into two sub-structures SS(1) with 11 internal DOFs and SS(2) with 8 internal DOFs, and 1-DOF of junction the mass m/2. The starting equation 20-DOF will be condensed and will bring to the resolution of a 10-DOF equation, divided into 1 junction DOF, 5 modes free or fixed interfaces of SS(1) and 4 modes free or fixed interfaces for SS(2).

Figure 1.

Decomposed structure.

The following characteristics are considered:

• Masses: $m_1=m_2=\cdots =m_{20}=2\mathrm{kg}$ ;

• Linear stiffness: $k_1=k_3=\cdots =k_{41}=50\mathrm{N}/\mathrm{m}$ ;

• Nonlinear cubic stiffness: $k_2=k_4=\cdots =k_{40}=k_{42}=10\mathrm{N}/{\mathrm{m}}^3$ .

The initial conditions are

$$\begin{array}{c}\left\{u_c\right\}=\{0,0,0,0,0,0.5,0,0,0,0\}\hfill \\ \left\{{\stackrel{·}{u}}_c\right\}=\{0,0,0,0,0,0,0,0,0,0\}\hfill \end{array}$$

To illustrate the steps of the previously presented method, one begins by writing the vibration of the overall system of equations and those subsystems.

The time responses of the entire system (without reduction) and the system decreased with a fixed interface and free interface (CMS) method are shown in Figures 2–4 which correspond, respectively, to temporal displacements of the mass (12) which corresponds to DOF junction, mass (1) of the sub-structure SS(1), and the mass (18) of the sub-structure SS(2). We can see that the different methods of modal synthesis provide very similar results. Other numerical results were obtained for different types of loading and modal decomposition. Adaptation to the assembled continuous systems is straightforward using the finite element method.

Figure 2.

The temporal displacement response for m(12).

Figure 3.

The temporal displacement response for m(1).

Figure 4.

The temporal displacement response for m(18).

In this study, it has been chosen to investigate the effects of uncertainties by considering mass uncertain parameters. The mass parameter is supposed to be a random variable and defined as follows: $m=m_0(1+{\sigma }_m{\vartheta }_m)$ where ${\sigma }_m$ is a zero mean value Gaussian random variable, $m_0=m_{ii=1,\dots ,20}=2\mathrm{kg}$ is the mean value, and ${\vartheta }_m=2\%$ is the standard deviation of this parameter. First, the mean and variance of the magnitude of displacement have been computed by the perturbation method for the entire system (without reduction). The obtained results are compared with those given by the direct Monte Carlo simulation with 700 samples. The obtained results are plotted in Figures 5–7 which correspond, respectively, to temporal displacements of the mass (12) which corresponds to DOF junction, mass (1) of the sub-structure SS(1), and the mass (18) of the sub-structure SS(2).

Figure 5.

(a) The mean of temporal response for m(12) and (b) the variance of temporal response for m(12).

Figure 6.

(a) The mean of temporal response for m(1) and (b) the variance of temporal response for m(1).

Figure 7.

(a) The mean of temporal response for m(18) and (b) the variance of temporal response for m(18).

Finally, we used the approach based on the coupling of the perturbation method with CMS condensation method. This approach allows reducing the size of the problem and the computational cost. The mean and variance of displacement have been shown in Figures 8–11. We can see that the different methods provide very similar results.

Figure 8.

The mean of temporal response for m(1), Monte Carlo simulation with 700 samples, and perturbation with whole structure and with free and fixed interface (CMS) methods, ${\vartheta }_m=2\%$ .

Figure 9.

The variance of temporal response for m(1), Monte Carlo simulation with 700 samples, and perturbation with whole structure and with free and fixed interface (CMS) methods, ${\vartheta }_m=2\%$ .

Figure 10.

The mean of temporal response for m(12), Monte Carlo simulation with 700 samples, and perturbation with whole structure and with free and fixed interface (CMS) methods, ${\vartheta }_m=2\%$ .

Figure 11.

The variance of temporal response for m(12), Monte Carlo simulation with 700 samples, and perturbation with whole structure and with free and fixed interface (CMS) methods, ${\vartheta }_m=2\%$ .

Conclusion

The main aim of this work is to provide the variability of the transient solution of a large and complex structure by considering geometric nonlinearities. We have achieved this by implementing an integrated approach, the coupling perturbation method, CMS reduction method, and temporal integration. The perturbation method was used to model the uncertain parameters by a development in Taylor series of second order. We developed the CMS method in the nonlinear case for reducing the finite element model. The implementation of the temporal integration by Newmark schema has allowed us to establish the variability of the solution for nonlinear reduced model with uncertain parameters. We could solve the problem of calculating the tangent matrix. The numerical tests show the accuracy of the results and minimization of cost calculation, thus validating this approach.