Nonlinear steady state vibrations of beams made of the fractional Zener material using an exponential version of the harmonic balance method

This paper presents the application of an exponential version of the harmonic balance method to the analysis of steady state vibration of geometrically nonlinear systems. A detailed description of the method and of the corresponding numerical procedure is provided. The von Karman theory is used to describe the effects of geometric nonlinearity. The material of the beams is modelled with the help of the Zener model using the fractional calculus. The problem is solved using an exponential version of the harmonic balance method. In the above-mentioned version, the complex calculus is used in contrast to the ordinary harmonic balance method, where the steady state vibrations are described with the help of the trigonometric functions. It significantly simplifies derivation of the amplitude equations. Moreover, the exponential version of the harmonic balance method is more elegant in comparison with the ordinary one. A detailed derivation of the amplitude equations is presented. The modified continuation method is proposed to solve the nonlinear amplitude equations and to determine the response curves. Moreover, the results of the exemplary calculation are presented and compared with known results in order to justify the efficiency and the correctness of the proposed approach.


Introduction
The harmonic balance method (HBM) is a very powerful tool, often used to solve problems of nonlinear vibrations. Using the HBM, the steady state vibration problems can be solved for strongly nonlinear systems [1][2][3]. Moreover, problems for which a multiharmonics solution was previously required could also be successfully solved [4,5]. The almost-period vibration is also analysed using HBM in [6][7][8]. The HBM is described in a number of papers [9][10][11][12] and monographs [13,14]. Due to the abundance of literature on the subject, only selected and most recent items are mentioned above.
In HBM, solutions to the equations of motion are sought (in almost all cases) which are in the form of trigonometric series (for one-dimensional discrete systems and in a case of forced vibrations) and can be written as where (t) is the vector of the sought solution of the equation of motion, ck and sk are unknown vectors (1) (t) = n ∑ k=0 ck cos k t + sk sin k t Abstract This paper presents the application of an exponential version of the harmonic balance method to the analysis of steady state vibration of geometrically nonlinear systems. A detailed description of the method and of the corresponding numerical procedure is provided. The von Karman theory is used to describe the effects of geometric nonlinearity. The material of the beams is modelled with the help of the Zener model using the fractional calculus. The problem is solved using an exponential version of the harmonic balance method. In the above-mentioned version, the complex calculus is used in contrast to the ordinary harmonic balance method, where the steady state vibrations are described with the help of the trigonometric functions. It significantly simplifies derivation of the amplitude equations. Moreover, the exponential version of the harmonic balance method is more elegant in comparison with the ordinary one. A detailed derivation of the amplitude equations is presented. The modified continuation method is proposed to solve the nonlinear amplitude equations and to determine the response curves. Moreover, the results of the exemplary calculation are presented and compared with known results in order to justify the efficiency and the correctness of the proposed approach.
By means of HBM it is possible to solve also undamped linear or nonlinear problems of free vibrations. In these cases is the natural frequency and only one trigonometric function (sine or cosine) is present in (1).
However, the second possibility is to use the complex representation of the Fourier series and to write a solution to the equation of motion in the form where k and k are unknown complex vectors, (•) is complex, conjugate and i = √ −1 is the imaginary unit. The exponential form of the solution of the equation of motion is very often used when linear dynamic systems are considered. This case is discussed in many monographs (see, for example [15]). In the context of mechanical nonlinear systems, the discussed exponential form of solution is almost non-existent. In open literature, the exponential version of HBM has been applied to some problems, particularly those related to nonlinear circuit analysis [16][17][18][19]. Only the application of the exponential version of HBM to a dynamic analysis of mechanical systems will be discussed. The most general and detailed description of this method is presented in [14]. The exponential version of HBM together with the continuation method was used in [11] to analyse periodic responses of rotor/stator dynamics. The discussed method was also utilized in papers [9,12] to discuss some problems of nonlinear dynamics. Using the exponential version of HBM, the steady state vibration is also analysed in [4]. In [4], the systems that feature distinct states (i.e. where the equations of motion involve piecewise defined functions) are considered. In paper [20], the application of exponential version of HBM to the analysis of vibration of bladed disks coupled by friction joints is presented. Complex nonlinear modal analysis for turbomachinery blades with friction and using the exponential version of HBM is presented in [21] by Laxalde and Thouverez. The application of HBM together with the continuation procedure is developed in [22] and applied to the analysis of the quasi-periodic solutions of the equation of motion.
In general, in order to obtain the unknown vectors ck and sk or k and k (for k = 1, 2, ..., n) the amplitude equations must be developed. This is usually done using the Galerkin method. The amplitude equations are nonlinear algebraic equations. It is common practice to determine the response curves. In this case, a family of solutions of amplitude equations are required for a set of excitation frequency values. The response curves are determined using the continuation method as a method for solving the amplitude equations. The abovementioned procedure is well known (and described, for example, in [8,23,24]) when the solution of the equation of motion is written in the form of trigonometric series (1). On the other hand, details of the exponential version of HBM and the solution to the equation of motion given in the form of Eq. (2) are not fully described and are not well known. In this case, as previously, the amplitude equations are nonlinear algebraic equations, though with complex coefficients this time. However, the derivation of the above-mentioned amplitude equations is somewhat simpler in comparison with derivation when the solution to the equation of motion is in the form of trigonometric series (1). It is especially important when the nonlinear parts of the equation of motion are analytic functions with continuous derivatives of the appropriate orders. The continuation method for the amplitude equation with complex coefficients has not yet been reported in the literature.
In this paper, the exponential version of HBM is presented together with a detailed derivation of amplitude equations for dynamical systems which nonlinear parts of the equation of motion are analytic functions. Moreover, the continuation procedure which can be used to solve the above-mentioned amplitude equations has been developed for a first time. Beams which execute steady state geometrically nonlinear vibration and made of a viscoelastic material are chosen as exemplary nonlinear dynamic systems. The Zener model with fractional derivatives is used to describe the mechanical properties of the beam's material. For a first time, the exponential version of HBM is applied here to analyse the steady state vibration of nonlinear systems (beams) where the fractional derivatives are used in the system description. The correctness and efficiency of the discussed approach is illustrated by the results of numerical calculation.
The paper is organized as follows. In Sect. 2, a viscoelastic (VE) beam is briefly described. The amplitude equations are derived in Sect. 3 and the continuation method which is applied to solve the equations of motion is presented in Sect. 4. The conditions of stability of the steady state solution are discussed in Sect. 5. The calculation results are presented and briefly discussed in Sect. 6 while the concluding remarks are presented in Sect. 7.

Description of nonlinear vibration of beams
The theory of VE beams, previously presented in [25], is briefly repeated here in order to make the presentation self-consistent. The beam with a total length L , area of cross-section A , mass m per unit length and the moment of inertia J is considered. The beam has its immovable end in the horizontal direction. Moreover, the horizontal inertia forces are neglected. The beam is loaded by the harmonically varying forces where is the excitation frequency, p c (x) and p s (x) are real "amplitudes" of excitation forces. Moreover, ) . The last form of description of excitation forces is more suitable for presented analysis.
The Euler-Bernoulli theory together with the von Karman theory are used to describe the influence of geometric nonlinearity. It is well known that the von Karman theory is valid when the amplitudes of vibrations are of the order of the beam's height h . A most important relations for considered beam are given in the Appendix 1.
After introducing Eq. (84) into Eq. (83), the following equation of motion is obtained which, however, cannot be written in terms of transversal displacement explicitly.
The beam is built of a VE material which is modeled with the help of the so-called fractional Zener model, which means that fractional derivatives are used in this model. As derived in [25], the constitutive equations written in terms of internal forces N(t) , M(x, t) and generalized strains (x, t) and (t) are: where is the relaxation time, E 0 is the relaxed elastic modulus, E ∞ -the non-relaxed elastic modulus. Moreover, the inequalities E ∞ > E 0 > 0 , > 0 must be fulfilled. The symbol D t (•) denotes the fractional derivative of the order α with respect to time t of the quantity in the parenthesis. In the description of rheological properties of VE materials the Riemann-Liouville definition of the fractional derivative is most frequently used [26,27]. However, di Paola et al. [28,29] have shown that in the context of rheology it is more logical to use the Caputo definition of the fractional derivative. Yet, it is known (see [27]) that for a system at rest at t = 0 or for systems that operate from t = −∞ , the Caputo fractional derivative is equivalent to the Riemann-Liouville derivative. In conclusion, when the above assumptions are fulfilled, the operator D t (•) could be understood in either way: as the Riemann-Liouville derivative or as the Caputo derivative. In order to be precise in the derivation the Riemann-Liouville definition is adopted in this paper. Thus, the fractional operator is given by: where is the gamma function. Moreover, it is important to note that if above definitions and assumptions are valid the fractional derivative of exponential function is D t (e i t ) = (i ) e i t (see [30]).
The application of fractional calculus to analysis of several mechanical problems is well-known. Some applications of fractional calculus in mechanics can be found in [31]. Moreover, the free, steady and transient vibration of beams, frame structures and plates where the fractional models are used to describe damping properties of systems materials are, among others, presented in [1,5,25,[32][33][34][35][36][37][38]. Some useful reviews papers which presented the application of fractional calculus in mechanics are [39][40][41].
Alternatively, the constitutive equations for the fractional Zener model can be written using the Boltzmann superposition principle. This approach leads to the equation of motion written in terms of displacements which is in the form of integral-differential equations. However, this possibility is not discussed here.
It is important to note that, in contrast to the purely elastic material and the Kelvin model of the VE material, the fractional derivatives of forces are on the left hand side of Eqs. (5) and (6). The consequence of this is that the equation of motion of beams cannot be explicitly written in terms of displacements (here, in terms of w(x, t) ). Instead, the virtual work equation is used which, for the considered beams, takes the following form: where the quantities with the symbol (•) are virtual ones.
The equations of amplitudes will be derived using the harmonic balance method. Because of this, the time averaged of the virtual work equation will be used and written as: where T = 2 ∕ is the period of excitation.

Derivation of amplitude equation using the exponential version of harmonic balance method
Only the one harmonic solution is assumed to describe the steady state vibration in this paper. As it is proved in many papers (see, for example, [13, , such solution accurately enough describes the dynamic behaviour of nonlinear systems when the systems are characterized by cubic nonlinearity and the secondary resonances are not present. The solution for steady state vibration is assumed in the following form: where the bar over the symbols means the symbol is complex, conjugate with an identical symbol without the bar. The functions for internal forces N(t) , M(x, t) and generalized strains (t) and (x, t) are chosen in such a way that the geometric relations (83) and (84) and the constitutive Eqs. (5) and (6) are exactly fulfilled. The mentioned solutions are assumed in the form: First of all, the amplitude equations for the continuous model of the beam are derived. By differentiating Eq. (9) with respect to the space variable and introducing the obtained result and the relation (13) into the nonlinear term of the motion Eq. (4), the following is obtained Next, relations (3), (9), (12) and (14) are substituted into Eq. (4) and the terms multiplied by T 1 (t) and T 1 (t) are separately equated to zero. Moreover, the terms multiplied by T 3 1 (t) and T 3 1 (t) are neglected. As the result, the following amplitude equations are obtained for the continuous beam's model: The above equations are particularly useful during the stability analysis of steady state vibration. Now, the amplitude equations for a discrete beam's model will be derived.
From the geometric relations (80) and (84), it is easy to obtain the following results: Relations between the coefficients of internal force functions and the generalized strain functions are obtained by substituting functions (10)-(13) into the constitutive Eqs. (5) and (6), and equating the coefficients on both sides of those equations standing near . Moreover, it must be taken into account that It is important to note this; the simple Leibniz rule is not valid for the fractional derivative (compare [43][44][45]). The final results are: Now the assumed solutions are introduced into the time averaged virtual work Eq. (8). The beam's virtual state is described as the following functions of time: where (•) means the virtual variation of (•) with respect to space and The subsequent part of the averaged virtual work equation will be analysed separately.
After introducing Eqs. (86)-(89) and (28) into the virtual work of excitation forces and performing integration with respect to time, the following results are obtained: In the course of calculation, the following integral appears: Proceeding in a similar way with the virtual work of inertia forces, the following is finally obtained: After substituting Eqs. (10) and (12) into the averaged virtual work of bending moment and integrating this part with respect to time, we obtain: Additionally, the above integral can be written in terms of "amplitudes" w 1 (x) and w 1 (x) when relations (19), (20) and (26) (13) and (27) are introduced into it and this term must be integrated with respect to time. During the calculation, integrals (31) are obtained. Additionally, the following is obtained It should be noted the terms in the averaged virtual work which give us the nonzero components of the amplitude equations are ones for which the sum of the exponents of functions appearing in integrals (31) and (36.1) is zero. The fact has greatly simplified derivation of amplitude equations.
Taking into account the above results of integration and remembering that (t) and N(t) depends on time only, the average virtual work of normal force can be written in the following form: Subsequently, using relations (21)-(23), (28) and (29), Eq. (37) can be rewritten in the form: where The discrete version of average virtual work equation is obtained using the finite element approximation for the "amplitude" functions w 1 (x) and w 1 (x) . The beam is divided into finite elements and in each element these functions are approximated using the well known cubic Hermite polynomials, i.e.
T are the vectors of nodal parameters and (x) is the matrix of shape functions, which can be easily found in many monographs. The virtual functions w 1 (x) and w 1 (x) are approximated similarly, i.e.
Using the well-known finite element procedure the following amplitude equations are obtained Detailed derivation of above equations is presented in the Appendix 3.
From the mathematical point of view, this is a set of nonlinear algebraic equations with respect to vectors and as the unknowns. However, in comparison with the usual cases, these equations have complex coefficients and complex, conjugate vectors are expected as the solution. This makes a significant difference in comparison with the ordinary HBM where the coefficients of the amplitude equations are real numbers. In many cases, the response curve is needed, which means that a set of solutions for different values of excitation forces are required. Therefore, the above amplitude equations are treated as equations with parameter and the main parameter is the excitation frequency .
Having the complex solution for the steady state vibration problem given by Eq. (9), this solution can be rewritten as a real function. This can be done using the Euler formula: and the final result is The response curve is determined using the continuation method. The continuation method, also termed the path following method, is frequently used to solve nonlinear equations with parameter, occurring in many problems encountered in modern mechanics. The static analysis of geometrically or/and physically nonlinear structures (see [46,47]) and the analysis of large-amplitude free and steady state vibrations [23,48,49] are some of the examples of such problems. A general description of the continuation method can be found, for example, in [50]. However, in all of the above-mentioned problems, a set of nonlinear algebraic equations with real coefficients must be solved. The expected solutions are also real vectors. The continuation method which will be used to solve the amplitude Eqs. (42) and (43) must be modified because they are a set of equations with complex coefficients. When the response curve must be determined, the vectors , and the excitation frequency is treated as unknown quantities.
The continuation method is an incremental-iterative one. The response curve is represented by a sequence of excitation frequencies and amplitude vectors, i.e., m , m , m , m = 1, 2, ... . For any incremental step, the vectors m , m and the frequency m of the preceding step is assumed to be known. The purpose of an incremental step is to find an increment of frequency Δ and increment of amplitude vectors Δ , Δ , which can be accumulated to yield Because the excitation frequency is treated as an unknown quantity, the number of the unknowns is greater than the number of equations. In order to ensure the uniqueness of the solution, some additional equations, called the constraint equation, have been added. However, according to the author's best knowledge, the continuation method was not applied previously to a set of nonlinear algebraic equations with complex coefficients. Owing to this difference in the types of coefficients the constraint equation must be defined differently.
The constraint equations in the following form are proposed: which is similar to the one proposed by Crisfield in [51], except that here it is adopted for a problem with complex vectors. It should be noted that the quantity s (the arch length) is a real number because the amplitude vectors are complex, conjugate.
The problem at hand is a nonlinear one and can be solved only by an iterative way. Let us assume that, after the iteration number i , we know an approximation of the solution denoted by (i) , (i) and (i) . The Newton method is used to determine the iterative change of the frequency increment and the amplitude increments and . According to the Newton method, the above mentioned increments are governed by the following equations: where the matrices qq , qq , qq , qq and vectors , are defined in the Appendix 4. The vectors of the amplitude changes and can be written as a sum of two components, i.e.: where 1 and 1 represent the influence of residual vectors , and the second terms (i.e., 2 , 2 ) are there due to . These vectors are obtained by solving Eqs. (50) with the appropriate right hand side of Eqs. (i.e. − , − when 1 and 1 are needed or − , − where 2 , 2 must be found).
The increment of frequency is calculated from the constraint Eq. (49). After substituting the total increments of and up to the (i + 1)th iteration given by into Eq. (49), the following quadratic equation with real coefficients and with respect to is obtained where The solutions to Eq. (53) will be denoted as 1 and 2 . Typically, they are real positive or negative numbers although complex numbers are also obtained in some cases. If the continuation method is applied to solve the nonlinear equation with real coefficients, the increments which give a positive angle between the previous and the current increments are chosen to avoid doubling back on the response curve. This condition could be written as Δ (i)T Δ (i+1) > 0 where Δ is (temporarily) the increment of amplitudes in the above-mentioned case. This condition cannot be directly adopted in the continuation method for the nonlinear equations with complex coefficients. The concept of angles in complex vector spaces presented in [52] is utilized instead. In this paper, the cosine of the angle c between two complex vectors and is defined as In the context of the problem at hand, the abovementioned angle between two successive increments of amplitudes is It should be noted that cos c is generally a complex number but the denominator in Eq. (54) is a real number.
At the end, the following criteria of choosing the right increment of frequency are proposed: However, some other cases must also be discussed. If, for both increments of , inequality (57) is fulfilled, the right increment of is the one which is close to lin = −c∕b . Moreover, if both increments, i.e., 1 and 2 give negative values of Re(Δ (i)T Δ (i+1) ) , then the current incremental step is restarted with the increment of the arch length Δs reduced to one half. The same procedure is followed in the case of complex solutions of Eq. (52). To prevent the number of iterations being too large, the maximum number of iterations is set. If the number of iterations exceeds the set maximum value, the incremental step is restarted according to the same procedure as before.
A new approximation of the solution to the amplitude equations is given by: The iterations are repeated until the inequalities are satisfied, where 1 , 2 and 3 denote the assumed accuracy of calculations.
At the end of the description of the continuation method, some remarks on the arch length increment and on how to start calculations are needed. A good initial approximation of the solution is needed at the beginning of each incremental step because the Newton method is only locally convergent. Normally, it is achieved when the solution obtained for a previous incremental step is chosen as the initial approximation of the solution in the current increment. But where a whole procedure is started, the initial approximation of the solution can be chosen as a solution to the corresponding linear problem calculated for the excitation frequency which is far from the resonance region.
The values of the arch length increment Δs could be a fraction of √ T , where denotes now the above-mentioned linear solution to the amplitude equations. A second possibility is to perform one or a few increments with the help of the Newton method and assuming that Δs = increment of amplitudes calculated for two successive values of excitation frequency. It is common practice to change Δs during the incremental process. After Crisfield [52], it is proposed to change Δs according to the following formula: where Δs p and I p are the arc-length and the iteration number in the previous incremental step, respectively, and I d is the desired number of iterations. In order to achieve enough points for good representation of the response curve, Δs is reduced (usually to one half) when the total increment of frequency Δ or the norm of total increments of amplitudes (i.e. √ Δ T Δ ) are greater than the assumed values.

Stability of steady state solution
The stability of the steady state solution is verified by the averaging method. The approach presented here is very similar to the one published in [30] but now the complex exponential functions are used to describe the steady state solution. The starting point of derivation is the equation of motion (4). Moreover, the small parameter 0 < ≤ 1 is artificially introduced, i.e., the term N(t)w ,xx (x, t) is written in Eq. (4) instead of N(t)w ,xx (x, t) . The steady state solution is assumed in the form of relationships (9)-(13) but w 1 (x, t) and where Equations (65) and (67) are now a set of first-order differential equations with respect to w 1 (x, t) and w 1 (x, t) and the solution to these equations is: According to the averaging method and taking into account the assumption concerning w 1 (x, t) and w 1 (x, t) , the average functions w 1 (x, t) and w 1 (x, t) could be used in the time range (t, t + T) instead of w 1 (x, t) and w 1 (x, t) . Moreover, it is possible to prove that ẇ 1 (x, t) =̇w 1 (x, t) and ẇ 1 (x, t) =̇w 1 (x, t) . The average equations are: It must be underlined that the R(x, t) function depends on w 1 (x, t) and w 1 (x, t) but, according to the average procedure, w 1 (x, t) and w 1 (x, t) are considered as being independent of time in the integration period. After substituting Eqs. (3), (9), (12), (13) and (68) into (70) and performing the necessary integration, the following equations are obtained: All the quantities with the wave which appear in Eqs. (71) above must be understood as average, slowly varying functions of time. Relationships between the averaged amplitudes of moments and normal forces and the averaged amplitudes of displacements are in the form of relationships (19)- (23).
Moreover, the averaged quantities are introduced in the place of their corresponding quantities without the wave.
When the steady state solution is considered, then ̇w 1 (x, t) = 0 and ̇w 1 (x, t) = 0 and Eqs. (71) has (for = 1 ) a form identical with the amplitude Eq. (15). Only the averaged quantities instead of non-averaged ones appear in Eqs. (71). It could be proved that, for = 1 , the discrete form of Eqs. (71) is identical with the matrix amplitude Eqs. (42) and (43). As in the previous relationships, the averaged quantities must be introduced in the place of the non-averaged ones. Now, it is possible to rewrite Eqs. (71) (for = 1 ) in the following matrix form: In order to examine the stability of the steady state solution (9) the in-time behaviour of a small perturbation of averaged amplitudes Δ̃ and Δ̃ will be analysed. As shown above, the unperturbed solution fulfils Eqs. (72) while the perturbed solution denoted as Fulfils a similar set of Eqs. (74) Now, Eqs. (72) is subtracted from Eqs. (74). Moreover, the nonlinear terms are expanded with respect to Δ̃ and Δ̃ into the Taylor series and only two terms of these series are retained. Finally, the following set of differential equations with respect to Δ̃ and Δ̃ is obtained: where qq , qq , qq and qq are the previously derived blocks of tangent matrices. These matrices depend on the amplitudes ̃ and ̃ of the currently examined steady state solution.
The It should be noted that a non-typical eigenvalue problem is obtained as the result. In contrast to an ordinary case, the first matrix in Eq. (77) is the matrix with complex coefficients while the second matrix is symmetric but with elements which are imaginary numbers. The above-mentioned eigenvalue problem needs a deeper analysis some time in the future but is out of scope here. However, some general remarks can be formulated: (i) If the real part of all eigenvalues are negative, then the examined steady state solution is stable; (ii) If at least one real part of the above-mentioned eigenvalues is positive, the steady state solution is unstable; (iii) At the limit between the stable and unstable solutions, the real part of at least one eigenvalue must be zero or both the real and the imaginary parts of eigenvalue are zero, i.e., = 0.

Discussion of results of exemplary calculation
The one-span VE beam subjected to the harmonically varying force acting in the middle of the beam's length is analysed. The beam is of the length L = 4.0 m and has the cross-section dimensions b = h = 0.4 m . The VE beam's material is described with the help of the fractional Zener model of which the parameters (taken from [30,53] First of all, convergence of the results for an increasing number of finite elements is verified. The response curves of a beam with fixed-fixed ends are shown in Fig. 1. The response in the vicinity of the first resonance region is presented. A variation of the amplitude of transverse displacement in the middle of the beam is shown. The response curve for 2 elements is represented by the solid line, for 4 elements-by the solid line with small rhombuses, for 6 elements-by the solid line with small circles, and for  Fig. 2. It should be noted that these parts look like coefficients multiplied by cosine and sine terms in a traditional form of the steady state solution. The influence of boundary conditions on the response curve in the first resonance region is illustrated in Fig. 3. It is easy to note that the influence of boundary conditions is significant; the largest nonlinearity effect is observed in the case of the simply supported beam while the smallest effect is seen for the fixed-fixed beam. These results are comparable In Fig. 4, the response curve for the simply supported elastic undamped beam (of which the elastic modulus is E 0 = 7.0 MPa ) is compared with the response curve for the VE beam. As expected, the response curve for the VE beam is almost perfectly limited by the response curve for the elastic beam. This limitation is not perfect because, in the fractional Zener model, the elastic properties of the material are present in the first term on the right side of Eqs. (5) and (6) but also in the second one. However, the influence of the second term is omitted from calculations for the response curve for the elastic beam.
Next, the response curves of a one-span simply supported-fixed beam are compared for various values of the fractional parameters . In Fig. 5, the plots of the response curves are presented for the fractional parameters = 1 (the solid line), = 0.8 (the solid line with circles) and for = 0.6 (the dashed line). Especially in the resonant range of the response curves, one can observe a significant influence of fractional damping on the progress of these curves. For = 1 , the maximum amplitudes of resonant vibration are much smaller than the maximum amplitudes of vibrations obtained when the order of the fractional derivative is = 0.6 . This indicates that reduction of the fractional parameters strongly increases resonant amplitudes. The effect is expected because, for a decreasing , the beam's material has lower damping properties. Differences in the progress of the response curves in the non-resonant range can also be observed although the deviations are not significant. The results presented in Fig. 5 are qualitatively consistent with the ones presented in [30] for the fixed-fixed beam.
The influence of relaxation time on the response curve for the simply supported beam is shown in Fig. 6. Three different values of the relaxation time are chosen, i.e. = 0.02 ms (the solid line); = 0.01 ms-the solid line with circles and = 0.001 ms-the dashed line. The maximal amplitude in the resonance region grows significantly for decreasing relaxation times.
The last parametric analysis shows the influence of the E ∞ ∕E 0 ratio on the maximal amplitudes of vibration in the vicinity of the first resonance region of the simply supported-fixed beam. The results of calculation are shown in Fig. 7, where the resonance curves of the simply supported-fixed beam are presented. The amplitude in that figure is the amplitude in the middle of the beam and the resonance curves are in the vicinity of the first resonance region. Four response curves are presented, i.e., the curve for E ∞ ∕E 0 = 1.43 shown as the solid line; the curve for E ∞ ∕E 0 = 2.0 shown as the solid line with circles; the curve for E ∞ ∕E 0 = 3.0 shown as the solid line with crosses, and the curve for E ∞ ∕E 0 = 4.0 shown as the dashed line. This ratio has a huge influence on the response. For the smallest E ∞ ∕E 0 ratio, the beam behaves as a strongly nonlinear system, whereas for the largest E ∞ ∕E 0 ratio, the beam can be treated as a nearly linear one. In the considered range of values of the E ∞ ∕E 0 ratio, the maximal resonance amplitude is a monotonically decreasing function of the E ∞ ∕E 0 ratio.

Concluding remarks
An exponential version of HBM, in which the steady state solution of the equation of motion is written using the complex exponential functions, is presented. The detailed derivation of the amplitude equations is described. The procedure for derivation of amplitude equations is much simpler than that of the ordinary version of HBM. Now, in contrast to the ordinary HBM, the amplitude equations are the nonlinear algebraic equations with complex coefficients. Such amplitude equations are treated as a set of equations with excitation frequency as the main parameter. For a first time, a new version of the continuation method (known also as the path-following method) is proposed for determining the response curves. Moreover, the proposed method is, for a first time, used for the analysis of systems for which the equation of motion contains fractional derivatives. The steady state, nonlinear vibration of beams is considered as an example of geometrically nonlinear systems. The von Karman theory is used to describe the behaviour of beams which are made of a viscoelastic material. The so-called fractional Zener model with the Riemann-Liouville derivatives is used to describe the mechanical properties of the material. To overcome the problems with the physical law involved, including the fractional derivatives of internal forces and generalized strain, the approach with time averaging carried out prior to the HBM and the space integration was proposed. The stability of the steady state solution obtained with the help of exponential functions is also presented. In the course of analysis of numerical results for several simple examples of vibrating beams, some conclusions were drawn. They are related to the level of nonlinearity observed in the response curves taking into account different values of the beam's parameters, such as the impact of boundary conditions, order of the fractional derivative, relaxation time, and ratio of material modulus. We believe that there are no apparent difficulties in the possible future application of the presented approach to the analysis of the steady-state vibrations of other systems, including laminated beams, plates or shells with viscoelastic layers.
Using the well-known finite element procedure, it is possible to rewrite the integral appearing in Eq. (33) in the following matrix form: (82) N ,x (x, t) = 0 (83) mẅ(x, t) = T ,x (x, t) + p(x, t) It is easy to find the discrete form of the remaining parts of the virtual work Eq. (30), which are: The global and matrices are well-known linear stiffness and mass matrices, respectively, and are defined as follows at the finite elemental level: 1 T The global v matrix is called the VE stiffness matrices and is defined (at the elemental level) as: Moreover, the and are complex, conjugate vectors of nodal excitation forces defined as follows at the elemental level: Because the virtual vectors and can be freely chosen from the average virtual work equation, the socalled amplitude Eqs. (32) and (33) are obtained.

Appendix 4
The matrices qq , qq , qq , qq and vectors , are defined as follow: Above block matrices of the tangent matrix are calculated from the recent approximation of the solution, i.e., (i) ,