In-plane vibration of curved beams subjected to moving loads using finite element method

In-plane vibration of curved beams caused by loads moving at constant velocity is studied by means of finite element method. Curved beam elements incorporating shear effects by means of first order shear deformation beam theory are used for this purpose. Damping effects are assumed to be negligible. The vibration response of the curved beam under loads moving at a constant velocity along the beam is found and compared with that obtained using ANSYS to validate the study.


Introduction
Curved structures are used extensively in engineering applications. Applications of practical interest where curved beams are subjected to moving loads are bridges and railways. Curved beams introduce more complexity in an analysis due to their coupled governing differential equations. In case of inplane vibration, the governing equations of motions are in the form of three coupled partial differential equations in three variables. When performing finite element analysis (FEA), it is a general practice to approximate curved structures using large number of straight beam elements [1]. Although this can yield good results for practical purposes, the inherent coupling that exists in a curved beam is not captured using straight beam elements. This problem can be eliminated by using curved beam elements. Curved beam elements can model a curved geometry exactly and if they are based on exact form of shape functions, only a single element is required for giving exact results in static analysis. In case of dynamic analysis, very accurate results are obtained using fewer elements as compared with straight beam elements [1]. The element matrices for curved beam elements can be directly assembled i.e. transformation matrix is not needed provided the radius of curvature of the structure is uniform and the formulation is done in local curvilinear instead of local Cartesian system [11].
Friedman and Kosmatka [1] and Litewka and Rakowski [3] both developed curved beam elements based on exact form of shape functions. To account for shear deformation, they used first order shear deformation beam theory. The developed elements do not exhibit shear and membrane locking effects. Yang et al. [7] used Galerkin method to derive analytic solutions for in-plane as well as out of plane vibration response of a curved beam subjected to moving loads. Li and Ren [2] also used Galerkin method to obtain analytic solutions. They considered moving three directional loads, warping resistance in torsion and inertia force. They used a rigorous set of governing differential equations incorporating the above effects. Both Yang et al. and Li and Ren neglected the effects of shear deformation which can be quite significant when considering thick beams. In this study, by means of finite element method (FEM) using curved beam elements that account for the effects of shear using first order shear deformation beam theory, the in-plane vibration response of a curved beam under the action of moving loads is found. The results from present work are compared with the results obtained from commercial FEA package ANSYS.

Governing Differential Equation
The curved beam is modelled using first order shear deformation beam theory. The three degrees of freedom (dof) of the curved beam are the axial displacement (u), radial displacement (v) and rotation of the cross-section ( ). The x coordinate is coincident with the axis passing through the centroid of the curved beam. The cross-section is assumed to be bi-symmetric i.e. the products of inertia are zero so that no coupling between bending and torsion exists. The kinetic energy (K) and strain energy (U) (which in this case is equal to the potential energy) of the curved beam is given by [1] (1a) where L is the length of the curved beam, R denotes the radius of curvature of the beam element, E is the elastic modulus, G is the shear modulus, A is the cross-sectional area, I is the second moment of the cross-sectional area and k is the shear correction coefficient that is used in first order shear deformation beam theory. The dot indicates differentiation with respect to time. The work done (W e ) by the distributed axial force (F), distributed radial force (J) and distributed moment (M) is given by

Finite Element Model
The curved beam element developed by Friedman and Kosmatka [1] is used for FEA in this study. The element has two nodes with three dof per node (u, v and ) as shown in figure 1. The local coordinate system of the curved beam element is denoted by and has its origin at the first node of the element as shown in figure 1. This element uses shape functions that exactly satisfy the homogenous form of static equilibrium equations. The homogenous form of static equilibrium equations are obtained by setting the acceleration and external force equal to zero in equations (2a-c).
Following Friedman and Kosmatka [1], a general solution to the differential equations (3a-c) can be assumed as follows (4a-c) Substituting equations (4a-c) in the equations (3a-c) results in three algebraic equations in the three coefficients (X, Y and Z) and the exponential parameter ( ).
(5a) (5b) (5c) For the existence of non-trivial solutions, the matrix formed by the coefficients of X, Y and Z in equation (5a-c) must be singular. This gives an equation to solve for the exponential parameter ( ) which after some simplification can be expressed as follows (6) From equation (6), the first solution of is zero which means a constant term is a solution of the differential equation. As is squared indicating multiplicity, x multiplied by a constant is also a solution. Corresponding to the expression in the parentheses, is the solution of which means trigonometric terms in is a solution of the differential equation. Once again as the expression in the parentheses is squared, x multiplied by the trigonometric terms is also a solution. A general solution must include all the solutions mentioned above. Hence, the variation of the displacements along the element is assumed as follows (7a) As the assumed variations of displacement are exact, substituting equations (7a-c) in the equations (3a-c) and requiring them to vanish identically, 12 constraint equations are obtained which leave only six constants independent. These six independent constants are obtained from the six boundary conditions. The constraint equations render some of the constants to be zero. Apart from the constants that are zero, any six of them can be chosen to be independent. However according to Friedman and Kosmatka [1], the constraint equations can be expressed in simpler form when the six chosen independent constants are as shown below in vector {d}. By using the constraint equations, all the constants can be expressed in terms of the six independent constants in a matrix form in the following manner The six independent constants can be related to the physical nodal dof using the boundary conditions. We substitute = 0 for node 1 and = for node 2 in equation (8) where denotes the length of the curved beam element and by using relation (9) (11b) Using the relations (8), (9),(11) and equations (1a-b), the kinetic energy and strain energy of the curved beam element can be written in a matrix form from which the mass matrix and stiffness matrix can be identified as follows The work done by the external moving load on the beam element can be written in a similar manner as the work done by the distributed load in the following way (13) where F, J and M represent the moving loads along the axial, radial and rotational dof respectively. The moving load on the element is represented using Dirac delta function as follows (14) where are the magnitudes and is shown in figure 2.

Moving loads
Consider n number of arbitrary loads each moving with different but constant velocities on a curved beam. Let the time at which the respective loads start acting on the beam be . We can then represent the above moving loads as represents the magnitude of the loads, represents their respective velocity and L denotes the length of the curved beam. The above representation of the moving load is with respect to the global coordinate system. To calculate the nodal force vector for the moving loads using equation (10), all the moving loads have to be represented in component form (i.e. radial, tangential and moment) in the local coordinate system of the element as done in equation (14). The relation between the two coordinate systems is shown in figure 2. When a load is moving over a beam, it will appear in different beam elements as time progresses. The moving loads will only contribute to the nodes of the elements on which they are present at any particular instant. Therefore, in the global force vector the nodes corresponding to the elements where the loads are present would be non-zero and remaining entries will be zero. An alternative to the above mentioned method particularly useful when very large number of loads act on the beam is to use the principle of superposition i.e. deflection of the beam due to multiple moving loads is the sum of deflections due to the individual moving loads.

Results and discussions
To validate the study, the time history of the midpoint displacements and rotation of the cross-section of the curved beam found using curved beam elements are compared with the corresponding values from ANSYS for both single and multiple (two) moving loads. Both ends of the beam are fixed as shown in figure 3. For single load case, the beam is subjected to a load moving at a constant velocity of 40m/s and a magnitude of 1kN as shown in figure 3. In case of two moving loads, both the loads move at 40 m/s with magnitude equal to 1 kN. The first load starts acting on the beam at time t = 0 and the second load starts acting from time t = 0.6 seconds. The material properties of the beam used in the study and its geometric parameters are as follows:-Elastic Modulus (E) = 210 GPa, Poisson's Ratio (µ) = 0.3, Cross-sectional out of plane thickness (b) = 0.01 m, Cross-sectional depth (d) = 0.6 m, Second moment of area (I) = 0.00018 m 4 , Radius of curvature of the centroidal axis (R) = 45.5366 m, Density (ρ) = 7400 kg/m 3 , Angle subtended by the beam (Φ) = 60 degree. FEA using curved beam elements is done by coding in MATLAB and the beam is meshed with 80 elements. The time response is obtained through direct time integration [4,5] using Newmark's constant average acceleration method with a time step size of 0.005 seconds. The curved beam in ANSYS is modeled as a structure in plane stress state. Element type Plane 183 which has quadratic variation in displacement is used. It can be used as a six node triangle element or as an eight node quadrilateral element as is done in the present study. The beam is meshed with 720 elements and 2413 nodes were formed. The beam thickness (b) is input as a real constant. The entire analysis of moving load is done through programming using ANSYS Parametric Design Language (APDL). It should be noted that rectangular coordinate system is used in ANSYS and the formulation of the curved beam element is done using curvilinear coordinate system. However, at the mid-point of the curved beam where the displacement is sought the Y-component of the rectangular system is same as the radial component of the curvilinear system (opposite direction) and the X-component is same as the tangential component. Figure 4 and figure 5 show the comparison of the midpoint radial and tangential displacement respectively with values obtained from ANSYS. From the comparisons, excellent agreement of the displacement values with those obtained from ANSYS is observed. To verify the correctness of rotation values obtained using curved beam elements, we compare the axial displacement of the point located on the top surface at the midpoint of the beam (Z) which is related to midpoint axial displacement and rotation with the values that are obtained from ANSYS. The axial displacement of the point located on the top surface at the midpoint of the beam (Z) is given by the following relation (17) where y is the distance from centroidal axis taken positive towards the centre of curvature.
(a) (b) Figure 4. Comparison of midpoint radial displacement (v) obtained using curved beam elements and ANSYS for single moving load (left) and two moving loads (right).
(a) (b) Figure 5. Comparison of midpoint axial displacement (u) obtained using curved beam elements and ANSYS for single moving load (left) and two moving loads (right). As the axial displacement (u) has already been verified and y takes a constant value (equal to -0.3 for the top surface) for a particular point, the correctness of Z depends on the correctness of . Figure 6 shows the comparison for the values of Z with the ones obtained from ANSYS. Once again there is an excellent agreement between the two compared values which confirms the correctness of . Thus, the verification of the midpoint displacements and rotations with the values obtained from ANSYS validates the present study. In the case of single moving load, the load leaves the beam at around 1.2 seconds and thereafter we have free vibration. As no damping was assumed, the magnitude of free vibration does not decay. In the case of two moving loads, we see that up to 0.6 seconds the response is the same in both figures (i.e. single load and multiple load case) as only one load is acting till 0.6 seconds and the magnitude of the load is the same for both cases. After 0.6 seconds, on closer observation of the displacements and rotation plots for both single and two moving loads, we can see that the resulting displacement/rotation in the latter case can be visualized to be the sum of the displacements /rotations due to the two individual loads acting alone i.e. the results can also be obtained by the principle of superposition. The first load exits the beam at around 1.2 seconds. After 1.2 seconds, the response is due to the second load and the residual vibration due to the first load. At around 1.8 seconds, the second load exits the beam and thereafter we have free vibrations of the beam subject to the initial conditions at the time when the second load leaves the beam. Again as damping is neglected the free vibration magnitude does not decay.

Conclusions
In this work, in-plane vibration response of a curved beam under the action of moving loads was studied by means of FEM using curved beam elements. Damping was not considered in this study. The time history of the mid-point of the curved beam was calculated and compared with the results obtained from commercial software package ANSYS to validate the study. The finite element model developed in this work is not only simple and easy to implement but also accurate and reliable and can be safely used at least as a first estimate by designers and researchers in the study of important practical moving load problems such as vehicle induced vibration of bridges.