CONTRIBUTION TO RANDOM VIBRATION NUMERICAL SIMULATION AND OPTIMISATION OF NONLINEAR MECHANICAL SYSTEMS

This study considered the solution of the stochastic vibration of nonlinear mechanical systems with Gaussian random excitations. It realised a short review of linearisations techniques in stochastic dynamics mainly with application in the area of truss finite element modelling. The presented method of statistical linearisation is applied to numerical testing. In the second part of the article, the sensitivity analysis of the first two stress statistical moments and structural weight minimising subjected to the random stress constrains presented by mean value and standard deviation was brought to the fore. Cross-sectional areas were used as optimising parameters.


INTRODUCTION
Modelling and analysis of dynamic properties of mechanical systems using linear computational tools are often the first approximation of the actual conduct of the investigated object. The obtained results must be considered with caution, and we must decide whether they are sufficient for us, or we will require more detail and complex study of the issue with the consequent formation of a more appropriate mechanical and mathematical model. If at least one member in the model is nonlinear, then the whole system acts as nonlinear. According to the effect of non-linear members on the dynamic properties, we distinguished systems with significant (strong) and a little significant (weak) nonlinearities. Of course, much more complex are tasks with significant nonlinearities, where the ambiguity of solution causes considerable problems.
The task difficulty may increase if the inputs to a nonlinear system are characterised as a random process or its parameters are random as well, and especially if the investigated processes are non-stationary and non-Gaussian [1,25]. Several different methodical approaches have been developed to solve such problems, applicable only to certain types of tasks. The conditions determining the applicability of individual methods are associated with the type of nonlinearity, its location, as well as the characteristics of random processes operating in the system [1,25].
The most common but simplest nonlinearity is the static transformation (non-inertial), characterised by the following equation (1) Equation (1) presents the functional dependency between the restoring force in a spring and the spring strain. A mechanical system with nonlinearity solved in a static sense is usually a rare phenomenon. In technical practice, we encounter special systems that do not allow neglecting the inertial and dissipative effects. This means that equation (1) in general is no longer a nonlinear algebraic equation, but a nonlinear differential equation. Kropáč in [14] divides the methods of solving such tasks into the methods of a local nature that well describes the behaviour of a system in a relatively small surrounding of a certain (working) point usually in equilibrium. There are known linearisation approaches that have experienced considerable success in various applications in the engineering practice. In contrast, there are methods of global nature that describe the behaviour of a system in a broader sense, that is, in the entire scope of variables. Solution of a specific task depends on the solvers and their ability to identify the nature of the problem and predict the solution.
As already mentioned, there are a number of independent approaches built on different principles that give acceptable results based on certain assumptions. The global approaches include analysis of nonlinear dynamic systems using the theory of Markov processes, which leads to solving the familiar equation by Fokker-Planck-Kolmogorov (FPK) [2,5]. In practice, we often encounter local approaches represented by the method of statistical linearisation in various modifications (Krylov, Bogoliubov, Caughey, Bolotin [2], Kazakov [12], Nigam [16], Roberts, Spanos [18,23], Elishakoff [7]), method of statistical quadratisation (Spanos, Donley [22]), functional method by Volterra and Wiener [14,24]. In addition to the above-mentioned methods, which have been the most discussed ones in the recent decade, other methods have been developed in the past, such as the asymptotic method by Krylov-Bogoliubov-Mitropolsky [2,12], suitably adapted methods of the small parameter, especially the perturbation version [16], harmonic linearisation [2] and mean values [2]. Thanks to computer technology also other methods have been given a green light, based on simple but timeconsuming and computationally intensive approaches. In particular, there are simulation methods in various modifications. Simulation approaches solve dynamic tasks directly in the time domain, which is demanding on computation time. If we add to the above the random nature of excitation and the need to obtain a complete picture of system behaviour for various excitation realisations, then we come to the Monte Carlo method [1,10,19,20]. Other authors have explored ways of increasing efficiency by combining the Monte Carlo method with the other aforementioned methods (for example, a combination of statistical linearisation and the Monte Carlo method [7]).

METHOD OF EQUIVALENT STATISTICAL LINEARISATION
Equivalent statistical linearisation (ESL) has been used for a relatively long period to deal with randomly excited nonlinear systems, especially in the frequency domain. It is an approximation method, in which solving a system of nonlinear differential equations is replaced by solving an equivalent linear system suitable to obtain the Fourier transform of the system. The method was first presented by Krylov and Bogoliubov [2,12], and further elaborated by Booton [1,2], Caughey [21], Kazakov [12], Iwan, Yang [22], Spanos [23], then evaluated by Roberts and Spanos in [18]. The development of ESL within nonlinear stochastic dynamics is also associated with the works of Foeter [23], Malhotra and Penzien [24], Iwan and Yang [21], Atalik and Utku [18], Iwan and Mason [24], Brücker and Lin [21], as well as many others. Especially, the 1960s to 1980s were a period of enormous increase in method applications.
We obtain the parameters of an equivalent linear system by satisfying a certain preselected criterion, under which we assess the conformity of the original and linearised model especially in the outputs that are characterised mainly by statistical moments. By their very nature, the presented techniques are a generalisation of deterministic linearisation methods by Krylov and Bogoliubov. The most frequently used criteria include:  criterion based on energy balance of an actual and equivalent model,  criterion of conformity of the corresponding mean values and response dispersions of the original and linearised nonlinear function with a random input (less suitable),  criterion of the minimum root-mean-square deviation of an actual and approximated (linearised, or squared) function.
It should be noted that, on one hand, the advantages of this method certainly include its simplicity and admissible estimate of the first two statistical moments, but on the other hand, the change in nature of the random variable in the output, for example, the Gaussian output does not correspond to the Gaussian input, as contemplated by the method's theoretical foundations. A significant shortcoming of ESL is the possible incompatibility in the spectral response of an actual and linearised system (for example, the difference in the characteristics of the function of power spectral density, etc.).
For a system with multiple degrees of freedom, we present a generalised solution methodology based on the procedure applied by Spanos and Donley in [20] to the method of equivalent statistical quadratisation. Let's consider a system of nonlinear differential equations in the form where M, B, K are the matrices of mass, damping and stiffness of the linear part of the system, x(t) is the vector of unknown displacements, f(t) is the vector of random excitation (7,8) where (*) is the symbol for a complex conjugate (matrix), Sww is the power spectral density of the function w(t) = 2() and is equal to 1, Hf () is the Fourier transform of the impulse function hf (t). The response centred component can be determined from the equation We replace this "original" nonlinear system with an "equivalent linear" system where ai are the equivalent linearisation matrices that we determine from the condition of minimum root-mean-square deviation of differences between the functional values of the actual and linearised function, that is Parameter e is thus the vector of residuals defined as follows By substitution of (12) to the condition (11) we obtain the system where the vector has dimensions 2n1, the symbol < > represents again the mean value operator, A is the matrix of linearisation coefficients a1 and a2 with dimension 2nn in the form 14) and the matrix of right-hand sides GN has dimension 2nn and the form (15) The system (13) represents n linear systems with 2n unknowns. The resulting equivalent system for centred component will be in the form (10).
In technical practice, we often encounter tasks where input to the system is the power spectral density (PSD) matrix of excitation. Then, we calculate the response according to known equations for linear systems, that is where the frequency response of a linearised system is and Sxx() will naturally be the PSD output (response) matrix, which determines the spectral properties of the system and by its integration we determine an estimate of second-order statistical moments (variance). The mean value is determined from (5). The actual implementation of the calculation can only be carried out iteratively, since the matrices a1 and a2 contain information about the first two statistical moments. The task leads to solving a system of 2n nonlinear equations, that is to (5) and where Dx is the matrix of variances (we are mainly interested in diagonal elements). It is usually challenging to perform the initial estimate of mx and Dx, therefore it is preferable to first carry out the calculation only with the linear part of the system, and then calculate the matrices a1 and a2 from the results obtained. These matrices then get more specific in each iterative step. The calculation is completed when the pre-defined convergence criterion is satisfied.

STATISTICAL LINEARISATION OF GEOMETRICALLY NONLINEAR ROD FINITE ELEMENTS
Our aim is to present possible ways of analysing geometrically nonlinear rod systems excited by random forces. First, for simplicity, we expressed the total Lagrange's formulation of geometrical nonlinearity of a plane rod (truss) according to [3,26]. We know that for relative elongation the following relation applies where ue = [u1, u2, u3, u4] is the unknown nodal displacement vector, the matrices BL and BNL have the form where a is the rod (truss) length. Internal forces in the element are expressed as follows Let V be the rod volume (V = const.),  is the normal stress in the element, Ke is the nonlinear rod stiffness matrix, fv = [fv1, fv2, fv3, fv4] T is the vector of internal forces, for which applies that fv3 = -fv1 and fv4 = -fv2. Then, the forces fv1 and fv2 will be where x = u3 -u1 and x = u4 -u2. Nonlinear relations (23) must be "linearised" in accordance with the statistical linearisation method, that is, create a linear equivalent model. We are building on the theoretical basis described in the previous subchapter. Let us resolve the vectors of displacement and internal forces into a unidirectional component and a centred random component, that is If we apply the criterion (11), we obtain a system of equations (13) that in our case is as follows We solve the system (25) numerically, which generalises the procedure also for application to other types of finite elements with the only difference that the dimension of matrices will change. In simpler cases, such as ours, (Kekv)e can be determined analytically, as presented in [3] by Cherng and Wen. They expressed the internal forces using the following equivalent relations with the ultimate notation of the equivalent relation between the internal forces vector fv and the displacement vector ue (28) The matrices Klin and KSL represent the linear and linearised component of the resulting stiffness matrix of plane rod element. We calculate the mean values fv1m and fv2m using the following equations (29) The linearisation coefficients will be as follows: information always tells us more about the impact or significance of each optimisation variable on the monitored mechanical quantities (displacements, stresses, accelerations, etc.). The sensitivity analysis process can be considered as a selection process for the subsequent calculation of significant optimisation variables [11,17]. It is the calculation of gradients by exact way or by apn proximation procedure. In the first case, the gradient or the gradient matrix (sensitivity matrix) is expressed exactly. In the second case, the relevant derivations are calculated using known numerical equations.
Nonlinearity brings considerable complexity to this process, especially in the case of analytical expression of the necessary derivation. Let us show how the problems in the case of stress sensitivity analysis of the truss element are mentioned.
We are searching a derivation of the average stress in the j-th element according to the i-th cross-sectional area (optimisation variable Xi) where the derivation of fv3mj = -fv1mj according to Xi is generally obtained by derivation of (29). The calculation of the stress dispersion derivation in the j-th element of i-th crosssectional area will be where the derivative of the variance Dfv3j according to Xi is generally obtained by deriving the matrix of dispersions of the internal forces of the j-th element Since (Kekv)j is a function of the displacement dispersion matrix (Du)j and the average values of the displacement vector (um)j, it is a very complex task. The centre of the presented analysis will be the calculation of the global vector derivation of the displacements average values um and the global matrix of displacements variance Du according to Xi. We obtain a system of nonlinear algebraic equations by derivation of state equations (5 and 34) The system of equations (39) can only be solved iteratively, which of course reduces the effectiveness of this approach. The advantage of exact procedures in sensitivity analysis, based on a smaller number of numerical analyses loses the timeliness, hence, numerical derivation is more appropriate.

NUMERICAL EXAMPLE
We will design the cross-sectional areas of the rod structure considering three optimisation groups of the same cross-section. The first group of elements consists of rods 8 to 14 with cross-sectional area A1, second rods 15 to 25 with area A2 and third rods 1 to 7 with area A3. The loading force F has a random character of frequency-limited "Gaussian white noise" with PSD: Sff = 15·10 6    Minor changes in stress distribution not only for average values but also for standard deviations can be determined from processed results. This was caused by the applied geometric nonlinearity as well as the optimisation process itself as documented in Tables 3  and 4. Consideration of nonlinear analysis computationally complicated the task, but results do not differ significantly from linear stochastic analysis. The reason may be that the geometric nonlinearity of oscillating rod structures does not appear significantly as in other cases.

CONCLUSION
Methods of nonlinear analysis and optimisation of mechanical systems are certainly a topical issue today. Many approaches based on simplifying assumptions of different importance have been suggested in the past. Some approaches, especially from the numerical point of view, have succeeded, some have lost importance. This also happened in stochastic nonlinear dynamics. Due to the strong hardware support, time-based Monte Carlo simulations [10] are now preferred, but the application of linearisation techniques remains a useful tool for frequency domain solutions. The approaches of statistical linearisation are an interesting alternative in connection with finite elements with geometric or material nonlinearity.
Structural optimisation is an essential part of creating and analysing virtual models of structures in construction, engineering or other industries. The results of the presented study confirm that the optimum design of nonlinear model parameters is specific compared to linear models. Nonlinearities can cause significant changes not only in the optimisation process itself [6], but also in the end results, overestimating or underestimating the first two statistical moments. In conclusion, the optimal design of a nonlinear mechanical system leads to the striking result that cannot be predicted in advance. Therefore, a basic linear analysis of simplified physical models alone cannot be relied upon in design. It is necessary to get as close as possible to reality and it is often non-linear.