JuliaFEM beam element implementation

This article describes implementations, theory, and usage of beam elements in JuliaFEM. JuliaFEM is an open-source finite element method solver written in the Julia programming language, see Frondelius et al. [1]. Julia itself is a programming language designed especially for numerical computing, including features like just-in-time code compilation and multiple dispatch of generic functions [2]. Basically, beam element can be thought of as a reduction of a full continuum model, which still preserves the most important characteristics of the underlying physical phenomena. In academia, beams have often been described in two dimensions to introduce the basic concepts of finite element method to students. However, since calculation models needed in practical simulation work are usually described in three dimensions, the physics of the beam elements must be also described in three dimensions. JuliaFEM also has more general model reduction techniques which are described in article written by Rapo et al. [3]. In this article, we introduce a practical example of a natural frequency analysis of a formula race car frame. The obtained results are compared to the corresponding results of the commercial software to validate the implementation. Similar comparison of JuliaFEM natural frequency solver with continuum elements is done by Rapo et al. [4], where the model also had contacts, see Aho [5]. A more complicated use case, where beam elements are used in the optimization loop to find the fuel pipe route that minimizes the given cost function, is presented by Rapo et al. [6]. JuliaFEM also has preliminary support for shells, see Niemi et al. [7]. For a comprehensive history of Wärtsilä’s industrial calculations, see


Introduction
This article describes implementations, theory, and usage of beam elements in JuliaFEM.
JuliaFEM is an open-source finite element method solver written in the Julia programming language, see Frondelius et al. [1]. Julia itself is a programming language designed especially for numerical computing, including features like just-in-time code compilation and multiple dispatch of generic functions [2].
Basically, beam element can be thought of as a reduction of a full continuum model, which still preserves the most important characteristics of the underlying physical phenomena. In academia, beams have often been described in two dimensions to introduce the basic concepts of finite element method to students. However, since calculation models needed in practical simulation work are usually described in three dimensions, the physics of the beam elements must be also described in three dimensions. JuliaFEM also has more general model reduction techniques which are described in article written by Rapo et al. [3].
In this article, we introduce a practical example of a natural frequency analysis of a formula race car frame. The obtained results are compared to the corresponding results of the commercial software to validate the implementation. Similar comparison of JuliaFEM natural frequency solver with continuum elements is done by Rapo et al. [4], where the model also had contacts, see Aho [5]. A more complicated use case, where beam elements are used in the optimization loop to find the fuel pipe route that minimizes the given cost function, is presented by Rapo et al. [6]. JuliaFEM also has preliminary support for shells, see Niemi et al. [7]. For a comprehensive history of Wärtsilä's industrial calculations, see Frondelius et al. [8]. Typical applications of beam elements in industrial simulation is to model pipes and bolts in locations where there are no high requirements for the accuracy of the results.
Theory of beams is given only briefly. For a more comprehensive description about the theory of beam, see e.e. Crisfield [9], Belytschko [10], Oñate [11], and Zienkiewich [12]. Beam implementations has been considered in this journal earlier by Paavola and Salolainen [13]. For the very first FEM-based beam element implementations, see articles written by Kapur [14], Archer [15], Thomas et al. [16], Reissner [17] and Cowper [18].  The difference between Euler-Bernoulli and Timoshenko beam theories: in Euler-Bernoulli beam theory, the planes normal to the midline are assumed to remain plane and normal, whereas in Timoshenko beam theory the planes normal to the midline are assumed to remain plane, but not necessarily normal.
Two types of theories are widely used in the development of beam elements. The most important kinematic assumption concerns the motion of the normals to the midline of the beam. In Euler-Bernoulli beam theory, the planes normal to the midline are assumed to remain plane and normal, whereas in Timoshenko beam theory the planes normal to the midline are assumed to remain plane, but not necessarily normal see figure 1. Timoshenko beam theory can be seen as an extension to Euler-Bernoulli beam theory, allowing the effect of transverse shear deformation. Timoshenko beam theory relaxes the normality assumption, and the motions of an Euler-Bernoulli beam are a subset of the motions allowed by Timoshenko beam theory [10]. To clarify this, let us study the equations of Euler-Bernoulli beam theory and Timoshenko beam theory in a quasistatic situation. The governing equations of Timoshenko beam are the following coupled system of ordinary differential equations: where L is the length of the beam, A is the cross section area, E is the elastic modulus, G is the shear modulus, I is the second moment of area, q (x) is distributed load and κ is called Timoshenko shear coefficient. Now, if the last term of equation (2) can be neglected, thus resulting and substituting this result back to equation (1) reveals the governing equations of Euler-Bernoulli beam theory: The main implication of the kinematic assumptions are that when using Euler-Bernoulli beam theory, transverse shear vanishes. For this reason, Euler-Bernoulli beam theory works well only in the case of slender beams, whereas Timoshenko beam theory, also known as shear beam theory, works well also in the case of thick beams.
Euler-Bernoulli beam has only single dependent variable whereas the Timoshenko beam has two dependent variables. Because second derivative of velocity appears in the expression for the rate-of-deformation, the velocity field must be at least C 1 for Euler-Bernoulli beam. That is, the shape functions of Euler-Bernoulli beam needs higher continuity requirements than Timoshenko beam, where it is enough to have C 0 continuity. According to Belytschko [10], the continuity requirement is the biggest disadvantage of Euler-Bernoulli and Kirchhoff-Love theories, since C 1 approximation is difficult to construct in multi-dimensions. In the study of plates and shells, theory corresponding to Euler-Bernoulli beam theory is Kircchoff-Love theory, having similar assumptions, and theory of Timoshenko can be replaced with the theory of Reissner-Mindlin, respectively.
In the derivation of beam element, we assume that beam is straight with the axis of definition in the x direction and the cross-section A in the yz plane. We also assume that the primary stress components are the normal stress σ x and shear stresses τ xy and τ xz acting on each cross-section. The rest of the stress components are either neglected or included as boundary loads to the beam.

Equilibrium equations
In the derivation of equilibrium equations, conservation of linear momentum and angular momentum must be considered. First we study the global form of a balance of linear momentum in the cross section A of the beam: which, given using index notation, is A ∂τ xy ∂x + ∂σ y ∂y + ∂τ zy ∂z The main assumption is that the the cross-section is either constant of varies slowly enough along the x direction so that we can use where S is the perimeter of the cross section. Thus we can write Rearranging terms allows us to write the equilibrium of linear momentum in a compact form ∂P ∂x where Here, P is axial force resultant and S y and S z are transverse shear force resultants, and q x , q y and q z are loads in x, y and z directions, respectively. In the same way than with the linear momentum, we can write the global form of the balance of angular momentum in the cross section of the beam as where x is some arbitrary point in cross section andθ is the angular velocity vector. We consider the right hand side first. By using vector identity A × (B × C) = (A · C) B − (A · B) C, we can write the rate-of-change of the angular momentum in a form which can be further expanded in Cartesian coordinate system to where the definition of the components of the inertia tensor are Writing the equation (18) using index notation and neglecting terms σ y , σ z and τ yz , (see also Zienkiewich and Taylor [12]), we finally end up to the following: Rearranging terms allows us to write the equilibrium of angular momentum in a compact form where

Kinematics
We assume that the cross-section of beam remains as plane in deformation process. Displacement field is expressed as where u, v, w are displacement of the axis defining the beam and θ x , θ y , θ z are small rotation angles about the coordinate axes. Strains are calculated from displacement field using small strain theory, According to the 3D beam theory, strain components ε y , ε z and ε yz vanish, ε y = ε z = ε yz = 0. The nonzero strain components are The main assumption in Euler-Bernoulli beam theory is that rotations θ y and θ z coincide with the slopes of the neutral axis: With this assumption, nonzero strain components are thus

Elastic constitutive relations
Linear constitution relation is assumed: which is expressed in component form as Constitutive matrix D can be deducted from the general 3D constitutive equation.

Weak form
To derive a weak form, we define so-called generalized local strain vectorε and its energetic conjugate, generalized stress vectorσ. For a generalized local strain vector, we aim to write the nonzero strain tensor components ε x , γ xy and γ xz in terms of first and second order partial derivatives of displacement field. Thus our relation looks like In the same way, we define a generalized stress vectorσ: By examining strain transformation matrix S 1 and stress transformation matrix S 2 , it can be seen that S 1 = S T 2 , which has importance in the following step. Internal virtual work can now be written as By defining a kinematic matrix B such thatε = Bu andσ =Dε, equation 46 can be written in a standard form A more detailed description about generalized quantities is given by Oñate [11].

Discretization and finite element approximation
For finite element approximation, translations and rotations must be approximated by appropriate basis functions. The simplest 3D Euler-Bernoulli beam has two nodes. Linear C 0 continuous interpolation can be used to approximate the axial displacement u 0 and for the twist rotation θ x . A cubic Hermite C 1 continuous approximation must be used for v c and w c . First we consider the C 0 continuous linear interpolation for axial displacement and twisting. Starting with a linear polynomial to find shape function N 1 (ξ), there exists coefficients α 0 , α 1 such that Substituting conditions (49) to equation (48), we get the following system of equations and solution for coefficients α is so the first linear shape function is Using the same approach for N 2 (ξ): there exists coefficients α 0 and α 1 , such that p 1 (−1) = 0 and p 1 (1) = 1, which results Therefore, linear interpolations for axial displacement and twisting are To construct C 1 continuous interpolation polynomials, similar approach is used. Because now the function and its first derivatives needs to be continuous, four conditions needs to be matched, thus the polynomial is cubic: and its derivative is It's important to notice, that the C 1 continuity must be fulfilled in physical coordinates and not in dimensionless coordinates. For that reason, derivatives must be taken with respect to potentially curvilinear coordinate of beam center line s, not ξ, using chain rule: dp 3 (ξ) ds = dp 3 (ξ) dξ dξ ds .
so apparently dp 3 (ξ) ds = dp 3 where is Jacobian determinant. Especially, in the case of isoparametric straight beam, where (e) is the length of the element e. Polynomial and its derivatives needs to be evaluated in the endpoints of the dimensionless coordinate system ξ 1 = −1, ξ 2 = 1, yielding Now, shape functions can be found in the same manner than with linear interpolation. Shape function M 1 (ξ) can be found by setting the first equation 64 to equal one, and rest of the equations to zero, that is, In the same way, proceeding each shape function at time, we find the functions: To calculate the derivatives of shape functions with respect to physical coordinates, chain rule is needed, see again equation 60. The first derivatives are and the second derivatives are The interpolations with bending shapes in v (ξ) and w (ξ) are then The displacement interpolation can be written conveniently using matrix notation as Generalized local strain iŝ Stiffness matrix isK Force vector isf Finally, matrices needs to be transformed to global coordinate systems with transformation where x 1 and x 2 are the coordinate vectors of the element end nodes. Unit vectors e 2 and e 3 are defined along the principal directions y and z at each node, respectively.

Usage example
JuliaFEM is used by running a JuliaFEM input file, which is a normal Julia script including all details of the analysis. In this usage example, an input file is presented and explained row by row. The example is a natural frequency calculation of a three-dimensional frame structure using beam elements. Only the frame geometry, including node and element sets, is imported from the Abaqus input file. All other necessary details are defined in the JuliaFEM input. The orientation of beam elements is a key factor when using beam elements. The orientations are defined by normal and tangent directions. The tangent is automatically calculated from coordinates, but the normal must be defined to the JuliaFEM input by the user. The example structure is a formula race car frame from Formula Student Oulu [19]. The Formula student Oulu race car is presented in Figure 2. The geometry of the example Figure 2: Picture of the Formula Student Oulu race car model is the same as that of the actual race car used by the team, but the cross-section and material values are simplified. All pipes are considered as steel pipes having diameter of 12 mm and thickness of 2.5 mm. The example frame model is presented in Figure 3. For pipes, the moment of inertia is the same in all directions. In such situations, functions and loops can be conveniently used in JuliaFEM inputs. When using packages in Julia, it must always be stated which packages will be used. FEMBeam is a package itself but it is included in JuliaFEM; thus, we only need to write using JuliaFEM. The first thing needed for this calculation is the mesh. The beam frame mesh in this case includes an Abaqus input file that we want to import to JuliaFEM. Mesh is imported using abaqus read mesh()-function from AbaqusReader, which is also a package included in JuliaFEM. Before reading the mesh, a file path to the location of the mesh file must be defined. All basic information about the mesh is imported to JuliaFEM, including node sets and element sets. The mesh is imported to variable mesh, which comprises a  The elements for the problem are created using create elements () -function. With inputs mesh and "FRAME", the element set called "FRAME" from the mesh is used to create beam elements. The material and cross-section values are defined using the update!() -function. In this case, all beam elements have the same values and they can be defined all at once. 5 beam_elements = create_elements(mesh, "FRAME") 6 update!(beam_elements, "youngs modulus", 210.0e3) 7 update!(beam_elements, "shear modulus", 80.77e3) 8 update!(beam_elements, "density", 7.800e-9) 9 update!(beam_elements, "cross-section area", 176.715) 10 update!(beam_elements, "moment of inertia 1", 11320.778) 11 update!(beam_elements, "moment of inertia 2", 11320.778) 12 update!(beam_elements, "polar moment of inertia", 22641.556) The update!() -function is also used to define normals for the beam elements. As previously mentioned, the moment of inertia is the same in every direction when the crosssection is circular. This means that the normals can be chosen freely as long as they form an orthogonal basis where one direction is the tangent of the beam. This for-loop finds a suitable normal direction vector for every beam element. In the term k, the function indmax() finds the index of the maximum element in a collection. In JuliaFEM, all different types of elements will be added to a problem that will be solved. The problem is created using Problem()-function; in this case, the problem type is Beam with the name "frame" and the third value given for the function is the number of degrees which in this case is six. The beam elements and boundary condition elements are added to the problem using the add elements!()-function. To calculate the natural frequencies for the beam frame, only one step is needed, which is modal analysis. The step is created using the Analysis ()-function. The number of calculated eigenvalues can be set by defining step . properties .nev. Problems are added to the step using add problems!()-function. Results can be written to a Xdmf file [20], which is a file format that can be processed with, for example, Paraview. The file is created using the Xdmf-function and in this example, it is created on the same folder with the Abaqus input file.

Results
The outputs of the println -functions presented above are as follows. Moreover, the eigenmodes were animated using Paraview and analyzed. With given boundary conditions, the eigenmodes of the usage example were just as predicted. The biggest displacement occurs on top of the main hoop of the race car. The first eigenmode is presented in Figure 4.
Natural frequencies were also calculated using a commercial program to obtain comparison results. The same input file was used with same material and cross-section values. The only difference in the input was the defining of the beam orientations, which was done by hand with the Abaqus CAE. Defining all beam orientations by hand is inconvenient and leads to some inaccuracy.
The five lowest natural frequencies from both JuliaFEM and commercial program results were compared. The differences in frequencies are under 1 Hz, which is caused by the difference in defining beam orientations. Moreover, the eigenmode animations from both JuliaFEM and the commercial program were compared. Both programs showed the same eigenmodes.

Conclusions
The implementation of beam elements in JuliaFEM was introduced with an usage example. The usage example model is a complex 3D beam frame with a high number of beam elements being oriented in almost every possible direction. The calculation results were compared to results from commercial FEM program. The results support each other, and it can therefore be stated that our implemented beam element works accurately.
One of JuliaFEM's biggest strengths is that JuliaFEM inputs are normal Julia scripts wherein the user can use functions and loops to get the task done. This strength is emphasized in calculations like the usage example presented herein where functions and for-loops can be used. Defining the beam element pipe cross-section orientations is just one good example of using scripting in JuliaFEM inputs. Using these types of functionalities may require some programming skills but the programming language used is simple and easy to learn even for first timers. One reason these type of usage examples are made is because studying these usage examples is an efficient way to learn the basics.