FACTA UNIVERSITATIS Series: Mechanical Engineering Vol. 19, No 4, 2021, pp. 781 - 804 https://doi.org/10.22190/FUME210324045W © 2021 by University of Niš, Serbia | Creative Commons License: CC BY-NC-ND Original scientific paper STUDY ON VIBRATION RESPONSE OF A NON-UNIFORM BEAM WITH NONLINEAR BOUNDARY CONDITION Peng Wang1, Nan Wu2, Haitao Luo3, Zhili Sun1 1School of Mechanical Engineering & Automation, Northeastern University, Shenyang, China 2School of Mechanical Engineering, University of Manitoba, Winnipeg, Canada 3Shenyang Institute of Automation, Chinese Academy of Sciences, Shenyang, China Abstract. Forced vibration of non-uniform beam with nonlinear boundary condition is studied in this paper by proposing an iterative model combining Adomian Decomposition Method and modal analysis. An exponentially tapered beam with a hardening nonlinearity spring boundary is simulated as a case study. The model accuracy is proved by comparing iteration results and analysis solutions with linear and weakly nonlinear boundary conditions. Sin-weep nonlinear frequency spectrum is then obtained by the proposed model. The influence of boundary nonlinearity on the vibration response of non-uniform beam is analyzed. And the effect of different excitation amplitudes on nonlinearity in the vibration response is studied. The mathematical model and numerical solutions proposed in this paper can be used to solve and analysis broad vibration problems on general non-uniform beams with different nonlinear boundary conditions under various excitations. Key words: Nonlinear boundary, Non-uniform beam, Iterative method, Adomian Decomposition Method, Duhamel integral, Vibration characteristics 1. INTRODUCTION In practical engineering applications, non-uniform beam structures, including functionally graded material (FGM) beam structures are widely used because they can optimize weight and change strength by changing cross-sectional area and material properties. Over the years, many experts and scholars have studied the dynamics of non-uniform beams [1-4], including vibration characteristic analysis (natural frequency and modal shape solution) and vibration utilization (energy harvester and stability analysis) [5-9]. At the same time, the dynamics of structures with nonlinear boundary (including multi segment linear boundary condition) is a hot topic that has been studied in recent years. In real life, spring [10,11], rubber bearing [12,13], Received March 24, 2021 / Accepted May 16, 2021 Corresponding author: Zhili Sun, Nan Wu Northeastern University, 110819, University of Manitoba, R3T 2N2 E-mail: zhlsun@mail.neu.edu.cn, nan.wu@umanitoba.ca mailto:zhlsun@mail.neu.edu.cn 782 P. WANG, N. WU, H. LUO, Z. SUN concrete and elastic material and foundation [14] and soil [15,16] all have nonlinear characteristics, which have been taken into account in the static analysis of structures, but the nonlinearity of boundary is usually ignored in the dynamic analysis of structures. For the vibration governing equation of beam structure, the boundary condition determines the result of solution. In order to obtain more accurate structural dynamic characteristics, it is necessary to consider the nonlinearity of boundary condition, so it is meaningful to study the dynamic characteristics of non-uniform beam structure with nonlinear boundary. At present, there are many methods to solve the dynamic problems of uniform beam structure with nonlinear boundary condition. Iterative method is one of the effective methods. Ma and Silva solved the fourth order differential vibration equation of uniform beam with nonlinear boundary by iterative method [17]. Sun and Wang studied the existence of monotone positive solutions for a class of elastic beam equations with nonlinear boundary conditions by monotone iterative method [18]. Dang and Huong reduced the nonlinear fourth-order problem to a series of second-order linear problems with linear boundary conditions, and proposed an iterative method for solving the fourth-order nonlinear equations of beam structures with nonlinear boundary conditions [19]. Liu and Li proposed a fast iterative method to transform ordinary differential equations into integral equations to solve nonlinear beam equations with nonlinear moment boundary conditions [20]. Alves et al. studied the existence of monotone positive solutions for a class of beam equations with nonlinear boundary conditions by monotone iterative method [21]. At present, the iterative method to solve the dynamic problem of uniform beam structure with nonlinear boundary condition is to solve the nonlinear equations about deflection. The deflection here has no practical physical significance, and it is more about solving a mathematical problem. In addition to the iterative method, the reproducing kernel method and the expansion method can be used to solve the nonlinear boundary problems of beam structures. In [22] and [23], the analytic approximate solutions of a class of fourth order differential equations with nonlinear boundary conditions are studied by using the iterative reproducing kernel method and reproducing kernel Hilbert space method, respectively. Geng and Cui obtained a series solution to solve the singular nonlinear second-order periodic boundary value problem in the reproducing kernel space [24]. Sedighi et al. redefined the preloading nonlinearity as the boundary condition of cantilever beam with a new exact equivalent function (EF) of preloading nonlinearity, and obtained the corresponding analytical solution by using the parameter-expansion method (PEM) [25]. Li and Zhang studied the existence and uniqueness of monotone positive solutions for a class of elastic beam equations with nonlinear boundary conditions based on a new fixed point theorem of generalized concave operators [26]. Sedighi et al. used the newly introduced equivalent function to model the preloaded nonlinear boundary conditions of the beam, and obtained the analytical solution of the nonlinear vibration equation of the beam by He’s parameter expanding method [27]. Wang et al. considered the nonlinear fourth-order two-point boundary value problem (BVP) of elastic beam equation, and studied the existence, nonexistence and uniqueness of convex monotone positive solution of elastic beam equation with parameter  by using the fixed point theorem of cone expansion [28]. In addition to the above methods, there are also a series of methods to solve the dynamic problems of uniform beam with nonlinear boundary. Song uses the theorem of infinitely many critical points to study the existence of infinitely many solutions to the boundary value problem of elastic beam deflection on a fourth-order nonlinear elastic foundation that depends on two real parameters [29]. Mao et al. studied the nonlinear response of flexible structures with nonlinear general support conditions by using the modal revision method [30]. Rahman et al. studied the forced nonlinear Study on Vibration Response of a Non-Uniform Beam with Nonlinear Boundary Condition 783 vibration of Euler Bernoulli beam on nonlinear elastic foundation by using the improved multilevel residual harmonic balance method [31]. Li and Xu proposed an accurate Fourier series method for vibration analysis of multi span beam systems under arbitrary boundary conditions [32]. Bai and Wang discussed the existence of positive solutions for a class of nonlinear fourth order beam equations by using fixed point theorem and degree theory [33]. Wu proposed an iterative numerical method to solve the exact forced vibration of a cracked beam by considering the multiple modes and bilinear characteristics of the cracked beam [34]. For the non-uniform beam structure with linear boundary, especially for the general non-uniform beam with both cross-section width and thickness changing along the beam length, Adomian Decomposition Method (ADM) is a very effective method to solve its natural frequencies and mode shape functions [5]. The specific solving process of ADM is to decompose the solution of the equation and express it in the form of infinite series sum [35]. It is not necessary to use linearization, perturbation, iteration, model simplification, difference method and finite element method to solve the vibration differential equation with nonlinear term. Keshmiri et al. has used ADM to do the research work on solving the nature (mode shape functions and nature frequencies) and energy harvesting of non-uniform beam with introducing Taylor series [5-7]. But it is hard to obtain the general vibration response solution of non-uniform beam with nonlinear boundary because the boundary varies with time during the vibration process for the nonlinear vibration system. Furthermore, the mode function and natural frequency also change with the boundary. Even though, there are also some literatures on the dynamic problem of non-uniform beam with nonlinear boundary condition. Based on the Hamiltonian principle, Lin derived the governing differential equations for the non-uniform time-varying elastic boundary conditions of the pre-twisted non-uniform beam with coupled bending vibration and solved them by the method of separation of variables [36]. Kuo and Lee et al. transformed the governing differential equations into a set of self-adjoint linear fourth-order ordinary differential equations with variable coefficients by using the perturbation method, and studied the static deflection of a general elastic end-constrained non-uniform beam on nonlinear elastic foundation bearing axial and transverse forces [37]. Lee et al. extended the Mindlin-Goodman method and used the exact solution of the general elastic constraint non-uniform beam given by Lee and Kuo to study the dynamic and static responses of the non-uniform beam with non-uniform elastic boundary conditions [38]. Tsiatas proposed a boundary integral equation method for solving nonlinear problems of non-uniform beams on nonlinear three parameter elastic foundation [39]. Jang proposed an analysis method of moderately large deflections, which effectively considered the geometric nonlinearity caused by the moderately large deflections and the non-uniformity of the beam, and successfully and completely solved the moderately large deflections of the infinitely large non-uniform beam on the base of nonlinear elasticity question [40]. Lohar et al. assumed that the beam is on an elastic foundation and bears uniformly distributed loads. Considering different boundary conditions, static and dynamic parts are used to solve the large-amplitude free vibration behavior of axially functionally graded beams with different tapers [41]. These methods are used to transform or approximate the vibration governing equations of non-uniform beam with nonlinear boundary. The transformation process is complex, and some transformed equations cannot be solved directly, so numerical approximation method have to be used to solve them anyway. From the literature review, different models and methods were proposed to study the vibration characteristics of uniform and non-uniform beam structure with nonlinear boundary, but to the best of authors’ knowledge, the solution without transforming the vibration 784 P. WANG, N. WU, H. LUO, Z. SUN governing equation to solve the vibration response of general non-uniform beam (linearly or nonlinearly tapered along both width and thickness of beam cross-section) with a nonlinear taper function describing the variation of both width and thickness along the length direction (including non-uniform cylinder beam) and nonlinear boundary has not been proposed. In this paper, the ADM method and an iteration process are introduced to solve, simulate and study the vibration response of non-uniform beam with nonlinear boundary condition. Under solid spring properties, the vibration response of the non-uniform beam can be solved by ADM and Duhamel integral. In the numerical example section, the influences of different excitation amplitudes and frequencies on boundary nonlinearity are studied. When the properties of the spring are determined, the proposed iterative method can be used to solve the vibration response of the non-uniform beam with nonlinear boundary. 2. THEORETICAL MODEL In this section, general mathematical model and numerical progress describing and solving the vibration of a non-uniform beam sitting on a non-linear boundary are presented. The nonlinear boundary condition is considered to be with multi-linear elastic properties, while ADM method is used to solve the natures of the non-uniformed beam with different linear elastic foundation, and the iteration numerical method is applied to solve the vibration response considering the linear boundary condition in each short time iteration step. 2.1. Vibration model of non-uniformed beam sitting on elastic foundation The equation of free vibration of a non-uniform beam without considering damping is given below, 2 2 2 2 2 2 ( , ) ( , ) ( ) ( ) 0 w x t w x t A x EI x t x x      + =      , (1) where x is the position variable along the beam length, ρ is the mass density, A(x) is the cross-sectional area, w(x, t) is the deflection function, t is time, E is the modulus of elasticity, I(x) is the second moment of area. The general boundary condition of the beam is defined as, 2 1 2 0 2 1 2 0 ( , ) ( , ) ( ) 0 ( , ) ( , ) + ( ) 0 L x T x w x t k w x t EI x x x w x t w x t k EI x x x = =     − =           =     , (2) 2 22 2 2 2 ( , ) ( ) ( , ) 0 ( , ) ( , ) ( ) 0 L x L T x L w x t EI x k w x t x x w x t w x t k EI x x x = =     + =           − =     , (3) Study on Vibration Response of a Non-Uniform Beam with Nonlinear Boundary Condition 785 where 1Lk and 2Lk are the tension spring constants on the left end, 1Tk and 2Tk are the torsional spring constants on the right end. In this study, a sample non-uniform beam structure with exponentially varied circular cross section sitting on a nonlinear tension spring foundation is analyzed and shown in Fig. 1. At the left fixed position, kL is a tension spring, which has a nonlinear property. Other non-uniform beam with general nonlinear boundary conditions can also be solved by the progress described below. With the multi-linear assumption of the non-linear spring stiffness, the vibration in each linear domain of kL can be solved by treating the system as a linear one. The whole nonlinear vibration response can be solved by iteration process with small iteration time step, while in each time step, the system is considered as linear with non-change spring stiffness constant. d Left end Free end d L 0 L L k Fig. 1 Main and right views of an exponentially increasing tapered beam with nonlinear boundary condition For the above sample structure, its boundary condition can be expressed as below, 2 2 0 0 ( , ) ( , ) ( ) 0 ( , ) 0 L x x w x t k w x t EI x x x w x t x = =     − =           =   , (4) 2 2 2 2 ( , ) [ ( ) ] 0 ( , ) ( ( ) ) =0 x L x L w x t EI x x w x t EI x x x = =   =            . (5) By using mode superposition method, w(x, t) in Eq. (1) can be decomposed into two parts, 1 ( , ) ( ) ( ) i i i w x t W x q t  = =  , (6) where Wi(x) is the i-th mode shape and qi(t) is the i-th corresponding generalized coordinate of the free vibration response or external force. Substituting Eq. (6) into Eq. (1), one ordinary differential equation is obtained as, 786 P. WANG, N. WU, H. LUO, Z. SUN 2 4 3 2 2 2 4 3 2 ( )( ) ( ) ( ) ( ) ( ) ( ) +2 + =0 ( ) ( ) ( ) i i i i i d I xdI x d W x d W x d W x A x W xdx dx I x I x EI xdx dx dx  − , (7) where ωi is the i-th natural frequency. In order to solve the above equation, ADM is applied [5]. The operator form of the Eq. (7) is rewritten as,   2 ( ) ( )'( ) ''( ) ( ) 2 '''( ) ''( ) ( ) ( ) ( ) i x i i i i A x W xI x I x L W x W x W x I x I x EI x  = − − + , (8) where Lx is the fourth order differential operator. Lx -1 is applied on the both sides of Eq. (8) at the same time, where Lx -1 is the fourth-order integral operator. 2 3 1 2 3 4 1 2 ( ) 2! 3! ( ) ( )'( ) ''( ) 2 '''( ) ''( ) ( ) ( ) ( ) i i x i i i x x W x C C x C C A x W xI x I x L W x W x I x I x EI x   − = + + +   − + −    , (9) where C1 to C4 are constants that can be determined by boundary conditions, the detailed progress defining C1-C4 can be found in [5]. Wi(x) is written in series form, 0 ( )= ( ) ii k k W x W x  =  , (10) where k is the number of terms in series form. The larger the value of k is, the more accurate the solution is. A precise solution is often obtained with very small values of k [35]. Substituting Eq. (10) into Eq. (9), we have 2 3 1 2 3 4 0 1 2 0 0 0 ( ) 2! 3! ( ) ( ) '( ) ''( ) 2 '''( ) ''( ) ( ) ( ) ( ) i i i i k k k k x k k i k k x x W x C C x C C A x W x I x I x L W x W x I x I x EI x    =    − = = = = + + +      − + −           . (11) For each term in the series, we can have 2 3 0 1 2 3 4 ( ) 2! 3!i x x W x C C x C C= + + + , (12) 1 2 1 ( ) ( )'( ) ''( ) ( ) 2 '''( ) ''( ) 0 ( ) ( ) ( ) i i i i k k x k k i A x W xI x I x W x L W x W x k I x I x EI x   − +   = − + −     . (13) Substituting Eq. (12) and Eq. (13) into Eq. (11), the i-th mode shape function, Wi(x), can be obtained. With different boundary conditions and structural design, the mode shape function from ADM can be different and hence not given here. The natural frequency can be obtained by introducing the mode function into the boundary condition and solving the Study on Vibration Response of a Non-Uniform Beam with Nonlinear Boundary Condition 787 eigenvalue problem. The detailed solving process of the mode shape functions and nature frequencies of the non-uniform beam with general linear boundary conditions can be found in references [5-7] and is hence not provided here. After the natures (natural frequencies and modes shapes) of the structure vibration are solved by free vibration analysis and ADM, the forced vibration response can be solved using the modal analysis. The forced vibration governing equation of non-uniform beam structure considering damping under action of F(x, t), which is from a base motion, Y sin(  t), is given below, 2 2 2 3 2 2 2 2 2 ( , ) ( , ) ( , ) ( ) ( ) + ( ) ( ) ( , ) =- ( ) sin( ) w x t w x t w x t A x EI x C x I x F x t t x x x t A x Y t          + =        , (14) where C(x) is the strain rate damping coefficient. Y is the amplitude of the base displacement and  is the angular frequency of the base vibration. Substituting Eq. (6) into Eq. (14), the following equation can be obtained, 2 2 1 2 22 2 2 2 2 1 1 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) + ( ) ( ) - ( ) sin( ). i i i i i i i i i d q t A x W x dt d W x dq t d W xd EI x q t C x I x A x Y t dtdx dx dx  =   = = +   =           (15) Based on the understanding from the free vibration governing equation, for the i-th mode of free vibration, we have 22 2 2 2 ( ) ( ) ( ) ( )i i i d W xd EI x A x W x dx dx     =    , (16) With the understanding from the relationship given in the above equation, both sides of Eq. (15) multiplied by Wj(x) (i=j) and integrated from 0 to L in space domain leads to, 2 20 1 2 22 2 2 20 1 1 2 0 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) + ( ) ( ) ( ) - sin( ) ( ) ( ) L i i j i L i i i i j i i L j d q t A x W x W x dx dt d W x dq t d W xd EI x q t C x I x W x dx dtdx dx dx Y t A x W x dx     =   = =   +     =     . (17) Since in general, the damping function does not possess the orthogonality property, it is assumed that the structural damping is in the form of C(x)=αE where α is a constant. According to the orthogonality of normal vibration modes, Eq. (18) can be obtained, while αωi 2=2ξiωi is defined, 22 2 0 2 2 0 - sin( ) ( ) ( )( ) ( ) 2 ( ) ( ) ( ) ( ) L ii i i i i i i L i Y t A x W x dxd q t dq t q t F t dtdt A x W x dx       + + = =   , (18) where ξi is the modal damping ratio of the corresponding i-th order natural mode. 788 P. WANG, N. WU, H. LUO, Z. SUN Eq. (18) can be solved by Duhamel integral, and the final time domain solution for i-th mode is, ( ) 0 2 0 5 2 2 2 3 2 2 2 2 2 0 2 2 2 2 2 2 2 2 1 ( ) ( ) sin( ( )) ( ) ( ) = ( (2 2 ) ( ) ) ( ) ( ) ( (( ) sin( ) 2 cos( )) (( ) si i i i i i i i i i i i i i i i t t i i d d L i L d i i d i i d i t i i d d d i i d d i i d q t F e t d Y A x W x dx A x W x dx e t t − − − =   − −  + − + + − − − + − + − +                                       n( ) 2 cos( ))) i i t t−    , (19) where id  is the damped frequency corresponding to the i-th vibration mode, 2= 1- id i i    , Fi() is the force coefficient corresponding to the i-th vibration mode. According to Eq. (12), Eq. (13) and Eq. (19), the response of the beam structure is, 1 ( ) 0 1 ( , ) ( ) ( ) ( ) ( ) sin( ( ))i i i i i i i t ti i d i d w x t W x q t W x F e t d          =  − − = = =   −    . (20) 2.2. Iteration process considering nonlinear boundary condition An iteration process is proposed in this section with multi-linear assumption for the nonlinear foundation of non-uniform beam. It is assumed that the elastic foundation in each short period iteration step is linear with constant stiffness constant, while the stiffness constant changes between different iteration steps. The premise of iterative method is that at t=0 s, the initial condition of vibration is known, and there is no force applied to the beam. At this time, the beam is considered to be at rest, and the deflection and velocity at any position of the beam are, ( , )=0 ( , ) =0 w x t dw x t dt      . (21) The step length of each iteration is Δt=tn+1-tn, 1≤n<∞, and n is the number of iteration steps. The initial condition for time t1 (t1=0 s) is, 1 1 = ( , )=0 ( , ) =0 t t w x t dw x t dt      . (22) Because some variables, like time, vibration natural frequencies, mode shape functions, in each iteration step are different in the process of iteration of the nonlinear system, they are defined by iteration step subscripts. The subscript 1 represents these variables in the first time period/iteration step in period t1-t2. The subscript 2 represents these variables in the second time period/iteration step in period t2-t3 and so on. The subscript n represents Study on Vibration Response of a Non-Uniform Beam with Nonlinear Boundary Condition 789 these variables in the n-th time period/iteration step in period tn-tn+1. Detailed specific presentation method of variables in the iteration section is shown in Table 1. i and j in Table 1 represents the i-th and j-th vibration modes, respectively. Table 1 Representation of variables in the iterative process Time period Variable (1≤i<∞, 1≤j<∞) t1≤t≤t2 1 ( , )w x t Total vibration response 1 ( ) i W x i-th mode shape function 1 ( ) i q t , 2 1 ( ) i t t dq t dt = i-th corresponding generalized coordinate of the external force and its derivatives 1i  , 1i d  i-th natural frequency and damping frequency 1 1 ( ) i F t + Force coefficient corresponding to the i-th vibration mode tn≤t≤tn+1 ( , ) n w x t Total vibration response ( ) in W x , ( ) jn W x i-th and j-th mode shape functions ( , ) nfree w x t Free vibration response in  , ni d  i-th natural frequency and damping frequency in A , in B Coefficients in the i-th corresponding generalized coordinate of free vibration ( ) in n F t + Force coefficient corresponding to the i-th vibration mode ( ) ni free n q t , ( ) ni n free t t dq t dt = i-th corresponding generalized coordinate of free vibration and its derivatives at time tn 1 ( ) in n q t + , 1 ( ) i n n t t dq t dt += i-th corresponding generalized coordinate of the external force and its derivatives at time tn+1 ( ) ni free q t , ( ) ni force q t i-th corresponding generalized coordinate of free vibration and external force during tn-tn+1 When the time is in t1-t2, there is no free vibration and the response of the forced vibration is, 1 1 1 1 1 1 1 1 1 1 1 - ( ) 1 1 10 1 ( , ) ( ) ( ) ( ) ( ) sin( ( )) i i i i i i i i i i d t t t t d w x t W x q t W x F t e t t d           = = − − − = = +    − −    ( 1 2t t t  ). (23) To simulate the nonlinearity of spring stiffness, kL is divided into n sections by different deflection intervals at the spring boundary location. It is considered that the value of kL corresponding to each deflection interval is the same, and kL, the mode shape functions and 790 P. WANG, N. WU, H. LUO, Z. SUN natural frequencies in different deflection intervals are different. By bringing each kL into the boundary condition, following the progress described in section 2.1, mode shape functions and natural frequencies can be obtained for the n-th time step (tn-tn+1) during iteration. The value of kL in the boundary condition of solving Wni(x) at time tn-tn+1 is the value of kL corresponding to the deflection at the spring boundary location, wn(x,t), x=0 for the sample structure shown in Fig. 1, at time tn. To clarify the iteration progress, we start the derivation from the second time step, t2-t3, while the first step vibration solutions is only from the Duhamel integral assuming the structure is at rest before excitation as described in Eq. (23). During the period, t2-t3, when the value of kL does not change from the one in the previous period, t1-t2, we have W1i(x) = W2i(x), 1i(x) = 2i(x), the free vibration response from the initial condition at t2 is, 2 2 2 22 2 2 1 ( , ) ( ) ( cos sin ) i i i i ii i t free d d i w x t W x e A t B t      − = =  + ( 2 3t t t  ), (24) where 2i A and 2i B are determined by 1 2( )i q t and 1 ( ) i dq t dt at t=t2. 2 2 2 22 2 2 2 1 2 ( cos sin ) ( ) i i i i ii i t d d e A t B t q t     − + = ; (25) 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 ( cos sin ) ( ) ( sin cos ) i i i ii i i i ii i ii i i t d d i d t d d i d t t t t e B dq t t t e A dt               − − = − − + = ; (26) 2i A and 2i B can be obtained by the above formula, 2 2 2 2 2 2 2 1 2 1 2 2 1 2 2 2 ( ) (( ( ) ) cos( ) ( ) sin( )) = i ii i i ii i i i i t i d d d t t d dq t e q t t q t t dt B         =  + + ; (27) 2 2 2 2 1 2 2 2 2 2 ( ) sin( ) = cos( ) i i i i i i i t d d q t e B t A t      − ; (28) During the period, t2-t3, the total vibration response is, 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 2 2 2 1 - ( ) 2 2 20 ( , ) ( ) ( )= ( )( ( )+ ( )) = ( )( ( cos sin ) 1 ( ) sin( ( )) ) i i i i i i i i i ii i i i i i i free force i i t d d i t t t t d d w x t W x q t W x q t q t W x e A t B t F t e t t d               = =  − = − − − =  + + +    − −     ( 2 3t t t  ). (29) where qfree2i (t) is i-th corresponding generalized coordinate of free vibration during t2-t3, qforce2i (t) is i-th corresponding generalized coordinate of the external force induced vibration during t2-t3. Study on Vibration Response of a Non-Uniform Beam with Nonlinear Boundary Condition 791 On the other hand, when the value of kL at time t2-t3 is different from that of t1-t2 with relative large boundary deflection, w(0, t), at t=t2 reaching to a different deflection interval, the corresponding vibration mode shape functions and generalized coordinates of free vibration at time t2-t3 also change. In order to calculate the generalized coordinates of free vibration at time t2-t3, the deflection function at time t2 need to be reassigned, 2 21 2 2 2 2 1 ( , ) ( , )= ( ) ( ) i i free free i w x t w x t W x q t  = =  (30) where w1(x,t2) is the deflection at time t2 from the previous iteration period, t1-t2.. wfree2(x,t2) is the response of the free vibration at time t2, 2 ( )i W x is the mode shape function of the structure starting at time t2, 2 1 ( ) ( ) i i W x W x , qfree2i (t2) is the i-th generalized coordinate of free vibration at time t2. Both sides of Eq. (30) are multiplied by 2 ( )jW x and integrated on 0-L giving 21 2 2 2 2 20 0 1 ( , ) ( ) ( ) ( ) ( ) j i ji L L free i w x t W x dx W x q t W x dx  = =   . (31) Considering the significant vibration natural frequencies and mode shapes changes due to the different boundary conditions with varying kL and the inevitable error in the mode shape functions derived by ADM method, when 2 2 ( ) i free q t (i=1,2,3) are solved, mode shape functions 2 ( )jW x (j=1,2,3) are introduced in Eq. (31), respectively. Through Eq. (31), if only the generalized coordinates of first three orders are taken consideration as an example, we can then obtain, 2 3 13 1 1 3 1 2 1 3 2 3 1 3 1 1 2 1 3 2 2 2 1 2 20 0 2 1 2 2 20 0 2 2 2 2 2 2 2 20 0 0 0 2 2 2 2 2 2 2 20 ( ) ( ( ) ( , ) ( ) ( ) ( , ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ) / ( ( ) ( ) ( ) ( ) ( ) ( ) L L free L L L L L L q t W x w x t dx W x dx W x w x t dx W x W x dx G W x W x dx W x W x dx G W x W x dx W x dx W x dx W x dx H W x W x dx W x W x dx = − + − −         1 3 2 3 1 0 0 0 2 2 2 2 2 2 20 0 0 ( ( ) ( ) ) ( ) ( ) ( ) ) L L L L L L L W x W x dx H W x W x dx W x dx− +        , (32) 2 22 3 2 2 ( ) ( ) free free q t G H q t= +  , (33) 2 1 2 1 21 2 2 1 33 1 2 2 1 2 2 2 20 0 2 2 20 2 20 ( ) ( ( ) ( , ) ( ) ( ) ( ) 1 ( ) ( ) ( ) ) ( ) L L free free L free L q t W x w x t q t W x W x dx q t W x W x dx W x dx = − −      , (34) where 792 P. WANG, N. WU, H. LUO, Z. SUN 1 2 1 2 1 2 1 1 2 2 2 2 1 2 2 2 2 1 20 0 0 0 2 2 2 2 2 2 20 0 0 ( ) ( ) ( , ) ( ) ( ) ( ) ( , ) ( ) ( ) ( ( ) ( ) ) L L L L L L L W x dx W x w x t dx W x W x dx W x w x t dx G W x dx W x dx W x W x dx − = −        , (35) 1 2 1 3 2 3 1 2 1 1 2 2 2 2 2 2 2 2 20 0 0 0 2 2 2 2 2 2 20 0 0 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ( ) ( ) ) L L L L L L L W x W x dx W x W x dx W x W x dx W x dx H W x dx W x dx W x W x dx − = −        , (36) 21 2 ( ) free q t , 22 2 ( ) free q t and 23 2 ( ) free q t are the first 3-th generalized coordinates of free vibration at time t2. We assume the free vibration response function during time t2-t3 to be, 2 2 2 2 1 2 2 2 2 2 1 ( , ) ( ) ( ) = ( ) ( cos sin ) i i i i i i i i i free free i t i w x t W x q t W x e A t B t      =  − = =   +   ( 2 3t t t  ), (37) where A2i and B2i here are determined by the initial conditions of the current iteration step at t=t2, 2 2 ( ) i free q t and 2 2 ( ) i free t t dq t dt = , which can be obtained from Eqs. (32-36). 2 2 2 2 22 2 2 2 2 ( cos sin ) ( ) i i i ii i i t d d free e A t B t q t     − + = ; (38) 2 2 2 2 2 22 2 2 2 2 2 2 2 2 2 2 2 2 2 ( cos sin ) ( ) ( sin cos ) i i i ii i i i ii i ii i i t d d i d freet d d i d t t t t e B dq t t t e A dt               − − = − − + = ; (39) A2i and B2i can hence be obtained by the above formula, 22 2 2 2 2 2 2 2 2 2 2 2 2 2 2 ( ) (( ( ) ) cos( ) ( ) sin( )) = i ii i i i i i i i i freet i free d d free d t t d q t e q t t q t t dt B         =  + + ; (40) 2 2 2 2 2 2 2 2 2 2 ( ) sin( ) = cos( ) i i ii i i i t free d d q t e B t A t      − ; (41) When kL changes from the previous iteration, during the second time step, t2-t3, the total vibration response is hence, Study on Vibration Response of a Non-Uniform Beam with Nonlinear Boundary Condition 793 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 2 2 2 1 - ( ) 2 2 20 ( , ) ( ) ( )= ( )( ( )+ ( )) = ( )( ( cos sin ) 1 ( ) sin( ( )) ) i i i i i i i i i ii i i i i i i free force i i t d d i t t t t d d w x t W x q t W x q t q t W x e A t B t F t e t t d               = =  − = − − − =  + + +    − −     ( 2 3t t t  ). (42) For general vibration solution during the iteration, when the time is tn-tn+1, if the value of kL is equal to the one in previous period tn-1-tn, the total vibration response is hence, 1 1 1 - ( ) 0 ( , ) ( ) ( )= ( )( ( ) ( )) = ( ) ( ( cos sin ) 1 ( ) sin( ( )) ) i i i n ni i i ni i i n i ni i n i n ni i ni ni n n n n free force i i t n n d n d i t t t t n n d n d w x t W x q t W x q t q t W x e A t B t F t e t t d               = =  − = − − − = +  + + +    − −     ( 1n nt t t +  ). (43) where ( ) ni free q t is i-th corresponding generalized coordinate of free vibration during tn-tn+1, ( ) ni force q t is i-th corresponding generalized coordinate of the external force during tn-tn+1. in A and in B are given as follow with known 1 ( )in n q t − from the previous iteration step, 1 1 1 ( ) (( ( ) ) cos( ) ( ) sin( )) = i d nn ii i i n n i ni i i n i ni t n i n n n d n d n n d n t t n d dq t e q t t q t t dt B         − − − =  + + ;(44) 1 ( ) sin( ) = cos( ) i n ni i i ni i ni t n n n d n n d n q t e B t A t     −  − ; (45) On the other hand, if the value of kL during the time period tn-tn+1 changes from that of tn-1-tn, the total vibration response is, 1 1 1 - ( ) 0 ( , ) ( ) ( )= ( )( ( ) ( )) = ( ) ( ( cos sin ) 1 ( ) sin( ( )) ) i i i n ni i i ni i i n i ni i n i n ni ni ni n n n n free force i i t n n d n d i t t t t n n d n d w x t W x q t W x q t q t W x e A t B t F t e t t d               = =  − = − − − = +  + + +    − −     ( 1n nt t t +  ). (46) Ani and Bni can be obtained by the above formula with the known vibration response, -1 ( , ) n n w x t , at t=tn from the previous step, 794 P. WANG, N. WU, H. LUO, Z. SUN ( ) (( ( ) ) cos( ) ( ) sin( )) = i d n nn ii i n n n n ni i i i i n i ni t free i n free n d n d free n d n t t n d dq t e q t t q t t dt B         =  + + ; (47) ( ) sin( ) = cos( ) i n ni n i ni i i ni t free n n d n n d n q t e B t A t      − ; (48) where the first 3-th generalized coordinates of free vibration at time tn, ( )nifree n q t (i=1,2,3), are given as below, 3 13 1 1 3 1 2 1 3 2 3 1 3 1 1 2 1 3 2 10 0 10 0 2 0 0 0 0 2 2 ( ) ( ( ) ( , ) ( ) ( ) ( , ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ) / ( ( ) ( ) ( ) ( ) ( ) ( n L L free n n n n n L L n n n n n L L L L n n n n n n n n n n n n n q t W x w x t dx W x dx W x w x t dx W x W x dx M W x W x dx W x W x dx M W x W x dx W x dx W x dx W x dx N W x W x dx W x W x − − = − + − −         1 3 2 3 1 0 0 0 0 2 2 0 0 0 ) ( ( ) ( ) ) ( ) ( ) ( ) ) L L L L L L L n n n n n dx W x W x dx N W x W x dx W x dx− +        , (49) 2 3 ( ) ( ) n nfree n free n q t M Nq t= + , (50) 1 1 21 2 1 33 1 10 0 0 2 0 ( ) ( ( ) ( , ) ( ) ( ) ( ) 1 ( ) ( ) ( ) ) ( ) n n n L L free n n n n free n n n L free n n n L n q t W x w x t q t W x W x dx q t W x W x dx W x dx − = − −      , (51) where we have, 1 2 1 2 1 2 1 1 2 2 1 10 0 0 0 2 2 2 0 0 0 ( ) ( ) ( , ) ( ) ( ) ( ) ( , ) ( ) ( ) ( ( ) ( ) ) L L L L n n n n n n n n n L L L n n n n W x dx W x w x t dx W x W x dx W x w x t dx M W x dx W x dx W x W x dx − − − = −        ,(52) 1 2 1 3 2 3 1 2 1 1 2 2 0 0 0 0 2 2 2 0 0 0 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ( ) ( ) ) L L L L n n n n n n n L L L n n n n W x W x dx W x W x dx W x W x dx W x dx N W x dx W x dx W x W x dx − = −        . (53) 3. NUMERICAL STUDIES, RESULTS AND DISCUSSION As shown in Fig. 1, a non-uniform cylindrical beam with positive exponential cross-section variation function is chosen for numerical case studies. The length of the beam is L=0.2 m and the diameter of the left end of the non-uniform beam is d=0.01 m. The diameter at x position along the length of the beam is ( ) n x d x d e  =  , with the taper ratio as n=2.0. The material of Study on Vibration Response of a Non-Uniform Beam with Nonlinear Boundary Condition 795 the beam is 6061 aluminum alloy with the elastic modulus as E =70×109 N/m2 and the density as ρ=2.7×103 kg/m3. The cross-section area and bending moment of inertia are 2 ( ) ( / 2) n x A x d e  =  and 4 ( ) 1/ 64 ( ) n x I x d e  =   , respectively. The sampling interval is 0.0001 s. The relationship between spring force and deflection is shown by a diagram in Fig. 2. In this paper, the continuous nonlinear boundary is treated as a multi- segment linear approximation. k1, k2, ..., in a diagram shown by Fig. 2 represent the kL values after the multi-linear segmentation. It is seen that the multi-linear assumption cannot present the exact nonlinear behavior/property of the reality nonlinear material with certain error leading to ‘slower’ hardening process with displacement increment. Considering certain number of linear segments, more obvious error can be noticed while the nonlinearity is more significant. However, such error can be further limited by introducing more refined multi-linear segments leading to more precise spring force considered during the calculation. In this work, we consider 10 segments for the case study to present the theory. The deflection intervals are determined according to the deflection range calculated by all kL values, and the left bound of deflection interval is defined to be smaller than the amplitude of steady state deflection. And it is more reasonable to define the interval according to the steady state amplitude to ensure that all kL values can be obtained during the vibration under harmonic excitation. When the deflection of the left end of the beam is greater than the maximum value in the interval of deflection, the elastic stiffness coefficient kL will change. There are many kinds of relationships between the tension (torsion) coefficient and deflection for different nonlinear springs, we just study one of them here. Deflection S p ri n g fo rc e Real deflection and spring force curve Deflection and spring force curve in calculation k1 k2 k3 k4 k5 Fig. 2 Relationship of deflection and spring force 3.1. Verification of mathematical model In this section, the accuracy of the proposed mathematical model is verified by using the iterative and non-iterative method to calculate the vibration response of non-uniform beam with linear boundary, and the iterative method to calculate the vibration response of beam with weakly nonlinear boundary defined by that kL changes within a small range 796 P. WANG, N. WU, H. LUO, Z. SUN (kL=20000-21800 N/m). The excitation angular frequency of  =672 rad/s and the base motion amplitude of Y=0.005 m are chosen for the verification. For the kL values with weak nonlinearity, the amplitudes of steady state deflections and their corresponding sectional intervals are shown in Table. 2. When values of kL are different, the mode functions and natural frequencies of the beam are different, and the amplitudes of steady state deflections at left end are also different. The amplitudes of steady state deflections at left end shown in Table. 2 will be larger closer to resonance. Since the variation range of kL is very small, it is considered that the values of kL change linearly with deflections variation at the left end of the beam. When the deflection at the left end is larger than the corresponding range of sectional intervals, the value of kL will change. Vibrations of the non-uniform beam at its free end in time domain are shown in Fig. 3. Table 2 kL values, amplitudes of steady state deflection and corresponding sectional intervals kL values (N/m) Amplitudes of steady state deflections at left end (m) Sectional intervals of deflection at left end (m) 20000 0.008154417514000 0-0.000855 20200 0.008205248588000 0.000855-0.00171 20400 0.008256681170000 0.00171-0.002565 20600 0.008308729902000 0.002565-0.00343 20800 0.008361411807000 0.00343-0.004275 21000 0.008414724816000 0.004275-0.00513 21200 0.008468682431000 0.00513-0.005985 21400 0.008523295153000 0.005985-0.00684 21600 0.008578582863000 0.0684-0.007695 21800 0.008634553316000 0.007695- Fig. 3 Deflections of non-uniform beam solved by iterative and non-iterative method in time domain Study on Vibration Response of a Non-Uniform Beam with Nonlinear Boundary Condition 797 It can be seen from Fig. 3 that the results of vibration response of non-uniform beam with linear boundary obtained by the iterative and non-iterative method are consistent with each other which preliminarily proves the accuracy of the iterative process. The difference between the vibration response of the beam with weakly nonlinear boundary and that with linear boundary solved by the iterative method is very small, which further proves the correct consideration of the nonlinearity with the iterative method. 3.2. Sin-sweep curves The vibration of the non-uniform beam with nonlinear boundary condition under different excitation frequencies are studied in this section. For the varying frequency of the base-motion excitation from 200 to 850 rad/s, the magnitude of the base-motion induced inertia force is kept equal. The kL values, the amplitudes of steady state deflections and their corresponding sectional intervals are the same as the values given in Table. 3. The value of kL is determined by the deflections at the left end and varies nonlinearly. It only represents one of the cases that the spring coefficient changes with the deflection. Maximum steady state vibration deflection at free end of non-uniform beam under a harmonic sweep test with different excitation frequencies are shown in Fig. 4. Table 3 kL values, amplitudes of steady state deflections and corresponding sectional intervals kL values (N/m) Amplitudes of steady state deflections at left end (m) Sectional intervals of deflection at left end (m) 20000 0.008154417514000 0-0.002 22759 0.008912748913000 0.002-0.004 31594 0.012528617870000 0.004-0.006 41120 0.020219170700000 0.006-0.008 50863 0.026634596390000 0.008-0.01 60702 0.017120586210000 0.01-0.012 70591 0.010843838260000 0.012-0.014 80511 0.007735876899000 0.014-0.016 90450 0.005971986427000 0.016-0.018 100402 0.004850948731000 0.018- It can be seen from Fig. 4 that for a beam with hardening nonlinear spring boundary (kL value changes following the deflection interval given in Table 3) under the increase of excitation frequency, the maximum steady state of the deflection will continue to increase and then suddenly drop. As a common nonlinear vibration phenomenon, this is because the natural frequencies of the beam with nonlinear boundary is varying during the vibration at different excitation frequencies and vibration amplitudes. Through the sin-sweep frequency progress, the deflection will become larger when the excitation frequency is close or equal to the higher natural frequency of the structure with hardening nonlinear spring boundary. When the excitation frequency is passing the range of the natural frequencies, the deflection will decrease significantly. At the same time, with the larger the excitation amplitude (0.0005-0.01 m), the stronger the level of the boundary nonlinearity can be noticed leading to wider frequency range with high vibration amplitude in the spectrum as shown in Fig. 4 (200-550 rad/s-200-765 rad/s). 798 P. WANG, N. WU, H. LUO, Z. SUN Fig. 4 Steady state vibration amplitude at free end of non-uniform beam under different excitation frequencies 3.3. Influence of the degree of boundary nonlinearity on non-uniform beam vibration Based on the proposed iterative method, further studies on the vibration response of non-uniform beam under nonlinear boundary condition are carried on. In order to study the influence of boundary nonlinearity on the vibration characteristics of non-uniform beam, the vibration response at the free end of the beam in time and frequency domain with fixed spring stiffness kL of 20000 N/m (linear boundary condition), 20000-31594 N/m, 20000-60702 N/m and 20000-90450 N/m are simulated. The excitation angular frequency is  =672 rad/s and the base motion amplitude is Y=0.005 m. Under this certain excitation, the multi-linear kL values, the amplitudes of steady state deflections at the left end of the beam and their corresponding sectional intervals are shown in Table. 3. The relationship between deflection and kL in Table. 3 is used for different degrees of nonlinear boundary. For kL in the ranges of 20000-31594 N/m, 20000-60702 N/m and 20000-90450 N/m, when the deflection is larger than the maximum value of deflection intervals (0.006 m, 0.012 m and 0.018 m) during vibration, the kL keeps unchanged. The vibration response at free end of the beam is shown in Fig. 5 (normalization was done in Fig. 5 (b)). The normalization process is to divide the magnitude by the maximum magnitude value to make the magnitude between 0-1. Study on Vibration Response of a Non-Uniform Beam with Nonlinear Boundary Condition 799 (a) (b) Fig. 5 Vibration response of non-uniform beam at free end with different degree of boundary nonlinearity ((a) time domain, (b) frequency domain) From Fig. 5, while the spring constant at the boundary is changing during vibration, compared with the linear boundary, with the increment of time, the deflection of the beam with nonlinear boundary will only approach a semi-steady state with variable amplitude in the time domain signal and clear wide bandwidth in frequency domain. This phenomenon is more obvious for stronger nonlinearity case with larger spring constant variation range (such as kL=20000-80511 N/m as given in Table 3). From the Fig. 5(b) it can be seen that for the 800 P. WANG, N. WU, H. LUO, Z. SUN beam with kL=20000 N/m, the frequency corresponding to the first peak and second peak are the natural frequency and excitation frequency. For the beam with kL=20000-31594 N/m, 20000-60702 N/m and 20000-80511 N/m, there are only one obvious peak. This is because the beam with linear boundary just has one natural frequency close to the excitation frequency. The natural frequencies of the beam with nonlinear boundary will change with the boundary conditions. The range of natural frequencies is relatively wide and includes excitation frequency, so from the Fig. 5(b), there will only be one peak at excitation frequency, and the peaks at natural frequencies are not very obvious. Fig. 6 Time periods of the system vibration staying at different kL values with different levels of boundary spring nonlinearity It is also interesting to see that when ranges of kL are 20000-60702 N/m, with the chosen base motion excitation constants,  =672 rad/s and Y=0.005 m, the beam experiences relatively larger vibration amplitude. But when the range of kL becomes 20000-31594 N/m or 20000-80511 N/m, the vibration amplitude decreases. This phenomenon is mainly due to the longer period of vibration close to resonance when the kL is between 20000 and 60702 N/m. From Table 3, it is noted that the beam steady state vibration amplitude reaches to relatively higher value, when the spring support stiffness, kL, is between 41120-60702 N/m leading to the natural frequency of the beam close to the excitation frequency,  =672 rad/s. Fig. 6 shows the times periods of the beam vibration staying at different spring support constants in different nonlinearity cases (in total 0.6 s). From Fig. 6, it can be seen that the beam experiences longer period of vibration at kL =60702 N/m closer to resonance, when kL is between 20000 and 60702 N/m compared with other two cases. Although the above values are just obtained from one specific case study, it can be concluded that if higher boundary of the hardening nonlinear spring stiffness leads to the resonance or close to resonance vibration, the vibration amplitude of the whole system can be larger, while the exact responses will differ with different system constants especially different nonlinear spring intervals. Study on Vibration Response of a Non-Uniform Beam with Nonlinear Boundary Condition 801 3.4. Influence of excitation amplitudes on boundary nonlinearity It is noticed that different excitation amplitudes may introduce different level of linearity (different nonlinearity spring stiffness variation range) in the vibration of the structures with nonlinear boundary. In this section, the excitation angular frequency is fixed as =672 rad/s. The kL values, the amplitudes of steady state deflections and their corresponding sectional intervals are shown in Table. 3. When the excitation amplitudes are 0.0005 m, 0.0025 m, 0.005 m and 0.01 m, respectively, the times periods of the beam vibration staying at different kL values (in total 1 s) are shown in Fig. 7. Fig. 7 Time periods of the system vibration staying at different kL values under different excitation amplitudes It can be seen from Fig. 7 that when the excitation amplitude is 0.0005 m, the beam vibration at the spring boundary is with extremely low amplitude (<0.002 m), and the spring stiffness does not change during the vibration with the linear boundary. With the increment of excitation amplitude, the range of the spring stiffness variation and level of the nonlinearity increases. When the amplitude becomes 0.01 m, the level of the boundary nonlinearity is the strongest covering all the kL variation range. At the same time, with different amplitudes, the vibration period of the beam staying at each kL is also different. Following the given the deflection intervals in Table 3, with different excitation amplitudes, the beam deflections at the spring boundary falling into each interval are different, and the vibration period on each kL is naturally not the same. When the spring properties are known, the vibration period of the beam staying at a certain kL can be controlled by changing the excitation amplitude, and then the vibration state of the beam can be adjusted. On the other hand, while the excitation frequency and amplitude are known, by adjusting the nonlinearity spring property (spring stiffness variation intervals), the vibration status can be controlled as well. 802 P. WANG, N. WU, H. LUO, Z. SUN 4. CONCLUSIONS In this paper, an iterative method is proposed to accurately solve the force vibration response of non-uniform beam with nonlinear boundary. Taking time node as iteration step, considering the change of natural frequencies and mode shape functions in the iteration process, the vibration response of beam in each short time period/iteration step is calculated by ADM and Duhamel integral considering both the transient and steady response. The influences of different excitation amplitudes and frequencies on boundary nonlinearity of non-uniform beam are studied. The results of numerical examples reveal the following conclusions: (1) Comparing iteration numerical results and analysis solutions with linear boundary condition, ADM method is found to be accurately combined with iteration progress solving non-uniform beam structural vibration with the nonlinear boundary condition, although inevitable minor error can be noticed in the calculated vibration mode shape functions from ADM leading to imperfect vibration modes’ orthogonality. (2) Through sin-sweep simulations, clear nonlinear spectrum with the hardening nonlinear boundary support can be noticed. With the increase of excitation frequency, the maximum steady state of the deflection will continue to increase and then suddenly drop because of non-resonance. Checking spectrums with different excitation amplitudes, the maximum steady state beam deflection will be larger and the frequency range with high vibration amplitude in the spectrum will be wider with the increase of the excitation amplitude. (3) Under the same base motion excitation, the beam vibration amplitude and the vibration period staying at a certain boundary stiffness range varies with different level of boundary nonlinearity. While the excitation frequency and amplitude are known, by adjusting the nonlinear spring property, the vibration status can be controlled. (4) Under the fixed nonlinear supported spring properties, the degree of nonlinearity reflected in the vibration response varies with different base excitation amplitudes. For vibration under harmonic excitation, the time staying at a certain range of boundary stiffness can be controlled by changing the base excitation amplitude. Acknowledgement: The study was funded by National Natural Science Foundation of China (No. 51775097). REFERENCES 1. Shabani, S., Cunedioglu, Y., 2020, Free vibration analysis of functionally graded beams with cracks, Journal of Applied and Computational Mechanics, 6(4), pp. 908-919. 2. Esmaeili, M., Tadi, B.Y., 2019, Vibration and buckling analysis of functionally graded flexoelectric smart beam, Journal of Applied and Computational Mechanics, 5(5), pp. 900-917. 3. Akbaş, Ş.D., 2019, Hygro-thermal nonlinear analysis of a functionally graded beam, Journal of Applied and Computational Mechanics, 5(2), pp. 477-485. 4. Chen, W.R., Chen, C.S., Chang, H., 2020, Thermal Buckling Analysis of Functionally Graded Euler-Bernoulli Beams with Temperature-dependent Properties, Journal of Applied and Computational Mechanics, 6(3), pp. 457-470. 5. Keshmiri, A., Wu, N., Wang, Q., 2018, Free vibration analysis of a nonlinearly tapered cone beam by Adomian decomposition method, International Journal of Structural Stability and Dynamics, 18(07), 1850101. 6. Keshmiri, A., Wu, N., Wang, Q., 2018, Vibration analysis of non-uniform tapered beams with nonlinear FGM properties, Journal of Mechanical Science and Technology, 32(11), pp. 5325-5337. Study on Vibration Response of a Non-Uniform Beam with Nonlinear Boundary Condition 803 7. Keshmiri, A., Wu, N., Wang, Q., 2018, A new nonlinearly tapered FGM piezoelectric energy harvester, Engineering Structures, 173, pp. 52-60. 8. Esmailzadeh, E., Ohadi, A.R., 2000, Vibration and stability analysis of non-uniform Timoshenko beams under axial and distributed tangential loads, Journal of sound and vibration, 236(3), pp. 443-456. 9. Keshmiri, A., Wu, N., 2019, Structural stability enhancement by nonlinear geometry design and piezoelectric layers, Journal of Vibration and Control, 25(3), pp. 695-710. 10. Sedighi, H.M., Shirazi, K.H., 2012, A new approach to analytical solution of cantilever beam vibration with nonlinear boundary condition, Journal of Computational and Nonlinear Dynamics, 7(3), 034502. 11. Wang, Y.R., Fang, Z.W., 2015, Vibrations in an elastic beam with nonlinear supports at both ends, Journal of Applied Mechanics and Technical Physics, 56(2), pp. 337-346. 12. Jankowski, R., 2003, Nonlinear Rate Dependent Model of High Damping Rubber Bearing, Bulletin of Earthquake Engineering, 1(3), pp. 397-403. 13. Ryan, Keri, L., Kelly, James, M., Chopra, Anil, K., 2005, Nonlinear Model for Lead–Rubber Bearings Including Axial-Load Effects, Journal of Engineering Mechanics, 131(12), pp. 1270-1278. 14. Rysaeva, L.K., Korznikova, E.A., Murzaev, R.T., Abdullina, D.U., Kudreyko, A.A., Baimova, J.A., Lisovenko, D.S., Dmitriev, S.V., 2020, Elastic damper based on the carbon nanotube bundle, Facta Universitatis-Series Mechanical Engineering, 18(1), pp. 1-12. 15. Trifunac, M.D., Todorovska, M.I., 1998, Nonlinear soil response as a natural passive isolation mechanism—the 1994 Northridge, California, earthquake, Soil Dynamics & Earthquake Engineering, 17(1), pp. 41-51. 16. Hartzell, S., Bonilla, L.F., Williams, R.A., 2004, Prediction of nonlinear soil effects, Bulletin of the Seismological Society of America, 94(5), pp. 1609-1629. 17. Ma, T.F., Da, S.J., 2004, Iterative solutions for a beam equation with nonlinear boundary conditions of third order, Applied Mathematics and Computation, 159(1), pp. 11-18. 18. Sun, J.P., Wang, X.Q., 2011, Monotone positive solutions for an elastic beam equation with nonlinear boundary conditions, Mathematical Problems in Engineering, 3, pp. 34-35. 19. Dang, Q.A., Huong, N.T., 2013, Iterative Method for Solving a Beam Equation with Nonlinear Boundary Conditions, Adv. Numerical Analysis, pp. 44-58. 20. Liu, C.S., Li, B.A., 2017, Fast New Algorithm for Solving a Nonlinear Beam Equation under Nonlinear Boundary Conditions, Ztschrift Für Naturforschung A, 72(5), pp. 397-400. 21. Alves, E., Ma, T.F., Maurício, L.P., 2009, Monotone positive solutions for a fourth order equation with nonlinear boundary conditions, Nonlinear Analysis Theory Methods & Applications, 71(9), pp. 3834-3841. 22. Geng, F., 2012, Iterative reproducing kernel method for a beam equation with third-order nonlinear boundary conditions, Mathematical Sciences, 6(1), pp. 1-4. 23. Geng, F., 2009, A new reproducing kernel Hilbert space method for solving nonlinear fourth-order boundary value problems, Applied Mathematics & Computation, 213(1), pp. 163-169. 24. Geng, F., Cui, M., 2007, Solving singular nonlinear second-order periodic boundary value problems in the reproducing kernel space, Applied Mathematics and Computation, 192(2), pp. 389-398. 25. Sedighi, H.M., Reza, A., Zare, J., 2011, Dynamic analysis of preload nonlinearity in nonlinear beam vibration, Journal of Vibroengineering, 13(4), pp. 778-787. 26. Li, S., Zhang, X., 2012, Existence and uniqueness of monotone positive solutions for an elastic beam equation with nonlinear boundary conditions, Computers & Mathematics with Applications, 63(9), pp. 1355-1360. 27. Sedighi, H.M., Shirazi, K.H., Reza, A., Zare, J., 2012, Accurate modeling of preload discontinuity in the analytical approach of the nonlinear free vibration of beams, Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 226(10), pp. 2474-2484. 28. Wang, W., Zheng, Y., Yang, H., Wang, J., 2014, Positive solutions for elastic beam equations with nonlinear boundary conditions and a parameter, Boundary Value Problems, 2014(1), pp. 1-17. 29. Song, Y., 2014, A nonlinear boundary value problem for fourth-order elastic beam equations, Boundary Value Problems, 2014(1), 191. 30. Mao, X.Y., Ding, H., Chen, L.Q., 2017, Vibration of flexible structures under nonlinear boundary conditions, Journal of Applied Mechanics, 84(11), 111006. 31. Rahman, M., Hasan, A.S., Yeasmin, I. A., 2019, Modified Multi-level Residue Harmonic Balance Method for Solving Nonlinear Vibration Problem of Beam Resting on Nonlinear Elastic Foundation, Journal of Applied and Computational Mechanics, 5(4), pp. 627-638. 32. Li, W.L., Xu, H., 2009, An exact fourier series method for the vibration analysis of multispan beam systems, Journal of computational and nonlinear dynamics, 4(2), pp. 710-733. 33. Bai, Z., Wang, H., 2002, On positive solutions of some nonlinear fourth-order beam equations, Journal of Mathematical Analysis and Applications, 270(2), pp. 357-368. 804 P. WANG, N. WU, H. LUO, Z. SUN 34. Wu, N., 2015, Study of forced vibration response of a beam with a breathing crack using iteration method, Journal of Mechanical Science and Technology, 29(7), pp. 2827-2835. 35. Adomian, G., 1988, A review of the decomposition method in applied mathematics, Journal of mathematical analysis and applications, 135(2), pp. 501-544. 36. Kuo, Y.H., Lee, S.Y., 1994, Deflection of nonuniform beams resting on a nonlinear elastic foundation, Computers & structures, 51(5), pp. 513-519. 37. Lin, S.M., 1998, Pretwisted nonuniform beams with time-dependent elastic boundary conditions, AIAA journal, 36(8), pp. 1516-1523. 38. Lee, S.Y., Wang, W.R., Chen, T.Y.F., 1998, A general approach on the mechanical analysis of nonuniform beams with nonhomogeneous elastic boundary conditions, Journal of Vibration & Acoustics, 120(1), 164. 39. Tsiatas, G.C., 2010, Nonlinear analysis of non-uniform beams on nonlinear elastic foundation, Acta Mechanica, 209(1-2), 141. 40. Jang, T. S., 1967, A general method for analyzing moderately large deflections of a non-uniform beam: an infinite Bernoulli–Euler–von Kármán beam on a nonlinear elastic foundation, Acta Mechanica, 225(7), pp. 1967-1984. 41. Lohar, H., Mitra, A., Sahoo, S., 2016, Geometric nonlinear free vibration of axially functionally graded non-uniform beams supported on elastic foundation, Curved and Layered Structures, 3(1), doi: 10.1515/cls-2016-0018.