Effects of meshing density of 1D structural members with non-uniform cross-section along the length on the calculation of eigenfrequencies

Density of division in finite element method does not affect only the accuracy of calculation, but also the necessary calculation time. This is directly influenced by the power of the used hardware, the efficiency of the algorithm used to assemble global stiffness and mass matrices and finally, by the method used to find eigenvalues of matrices for determination of eigenfrequencies. In engineering practice, when commercially available software is used, it is necessary to look for the optimum between the accuracy of the calculation and the length of the calculation. This paper deals with solution of eigenfrequencies of 1D elements with nonuniform cross section using Python 3.7.4 with libraries "numpy 1.18.1" for finding eigenvalues of matrices and "scipy 1.4.1" for finding solution for system of nonlinear equations.


Introduction
In finite element analysis, choosing mesh density is the finding optimum between accuracy of solution and complexity of model. In some of the simplest static analysis tasks, it is possible to effectively use the coarsest meshing. More complicated static analysis tasks require a finer meshing of the structure. The issue is discussed in more detail in the articles [1] and [2]. The influence of the meshing density of the structure on the results of modal analysis is discussed in the article [2]. Following chapters describe creation of a finite element method program in Python 3.7.4, analyse influence of dividing density on accuracy of calculation of the eigenfrequencies of beam with non-uniform cross-section. Program was designed to be able simply increasing number of finite elements in discretization of the simply supported beam using a single variable. To verify correctness of the algorithm of composing the global stiffness and mass matrices, the results were compared with an analytical solution of the calculation of the eigenfrequencies of the simple beam.

Programming of FEM
As it was mentioned in the introduction, for the purpose of analysis in the environment of Python, the FEM program was created. Subsequently, the "eigenvals" commands were called to obtain a vector of eigenvalues whose elements, after square rooting and dividing by τ represent the individual eigenfrequencies of the structure. To verify the functionality of algorithm which assembles global matrices of task, numerical calculations were compared with results of the analytical solution.

Principles of FEM
The first basic step of finite element analysis of structures is the discretization of the structure to finite elements. Since the task was to monitor the effect of division on the accuracy of the calculation. The next step is to approximate the deformation variables on each finite element separately. We start from the stiffness matrix of the finite element of the beam ( Fig. 1) with two degrees of freedom in both nodes. As it is shown in Fig. 1 in this case, chosen finite element has two nodes with three unknown deformations in each of them, thus stiffness matrices of elements (1) and their mass matrices (2) as well are 4x4 dimension [3] where: m i is the mass matrix of the i th element μ i is the mass per unit length of the i th element Next step is assembling global matrices of stiffness and mass from element matrices, respecting interconnections between individual members. Procedure of assembling global matrices is well known and it is identical for both, global stiffness matrix and global mass matrix. Scheme of the task is shown in Fig. 2 where the beam is divided into 4 finite elements.

Fig. 2. Scheme of the 4-element simple beam.
To obtain the eigenfrequencies of the structure, an eigenvalue function is called. Eigenvalue function in Python [5] is based on Linear Algebra Package, originally written in FORTRAN 77. In Python code it is necessary to import "numpy" library first. Then it is possible to call function for computing eigenvalues.

Eigenfrequencies of a simple beam
For purpose of verifying the functionality of the program, results of the numerical calculation method were compared with the results of the analytical solution. A simple beam, with uniform cross-section along its length and uniformly distributed mass, was solved according to the following scheme (Fig. 3).

Analytical solution
Based on [4], the analytical solution for a simple beam with a uniformly distributed mass is based on a partial differential equation describing the transverse damped beam oscillation with the uniformly distributed mass.

Numerical solution
The calculation of the eigenfrequencies of the simple beam was performed on the models with gradual increase in the number of finite elements of the given structure. Specifically, the beam was divided in range of finite elements from 4 to 1000. Results of numerical solution quickly converged to results of analytical solution, thus only results of division up to 100 finite elements are relevant.

Results comparison
The following chart (Fig. 4) shows the deviation of the numerical solution from the analytical solution for individual eigenmodes. The rapid convergence to the analytical solution in increasing the number of finite elements is evident. Coarse mesh in this relatively simple task is effective for obtaining results of first few eigenfrequencies. The graph shows a clear dependence -the higher the eigenmode shape, the higher the number of elements needed. Results of numerical solution quickly converge to results of analytical solution. Division of beam to more than 100 finite elements is practically noneffective, thus unnecessary.

Eigenfrequencies of non-uniform beams
Following subchapters deal with solutions of simply supported beams with variable crosssection height along its length. All beams have the same span of 12 m. Since exact solution of equation (3) is known only for prismatic simple beam, the numerical finite element method was used for the analysis of the beams. Again, the same finite element divisions were used as in Chapter 3.2. This time the results could not be compared with the calculated eigenfrequencies, thus values when beam was divided to 100 finite elements were used as a reference.
For beams with variable height of cross-section along its length, the function of height must be defined. The cross-sectional characteristics and uniform mass divided into nodes are calculated based on the average cross-sectional height between nodes of finite element. For each type of variable cross-sectional height, two beams with different steepness of the function defining height along its length were considered. The Young modulus of elasticity with value of 32.5 GPa was considered in all cases.

Variable height with the shape of a quadratic function
Variable cross-section height was defined by function (8)

Variable height with the shape of a cubic function
Variable cross-section height was defined by function (10) in case of beam 4a and by function (11) in case of beam 4b. Again, results are shown in the following charts (Fig. 7).  Fig. 7. Deviation of results depending on the division density of beams 4a and 4b.

Conclusion
As it was mentioned in the introduction, computational time is, among the others, affected by the complexity of algorithm which assembles global matrices. In the case of the thirddegree complexity of the algorithm, the computational time increases cubically with the number of finite elements. As it can be seen in the graphs in Chapter 4, increasing the dividing density of beams with variable cross-section is only effective to some extent. From these graphs it can also be seen that the steepness of the function describing the crosssectional height along the length of the beam influences the accuracy of the calculation. It can be stated that the division into 20 finite elements proved to be sufficiently precise and at the same time effective with respect to the necessary calculation time.