Implementing model reduction to the JuliaFEM platform

Summary. JuliaFEM is an open source (cid:28)nite element method solver written in the Julia language. This paper presents an implementation of two common model reduction methods: the Guyan reduction and the Craig-Bampton method. The goal was to implement these algorithms to the JuliaFEM platform and demonstrate that the code works correctly. This paper (cid:28)rst describes the JuliaFEM concept brie(cid:29)y after which it presents the theory of model reduction, and (cid:28)nally, it demonstrates the implemented functions in an example model. This paper includes instructions for using the implemented algorithms, and reference the code itself in GitHub. The reduced sti(cid:27)ness and mass matrices give the same results in both static and dynamic analyses as the original matrices, which proves that the code works correctly. The code runs smoothly on relatively large model of 12.6 million degrees of freedom. In future, damping could be included in the dynamic condensation now that it has been shown to work.


Introduction
JuliaFEM is an open-source nite element solver written in the Julia programming language [1].It enables exible simulation models, and it takes advantage of the scripting language interface which makes it easy to learn and embrace.It is a real Julia language programming environment where Structural Analyst or Researcher can combine FEM simulation with other analyses and work-ows.These features introduce an open source platform for testing new ideas and simulation models to the academic world in the nite element domain.[2,3] The JuliaFEM is installable meta-package, which means it is a collection of other Julia packages, which it installs as a dependency.ModelReduction.jl is one of the sub packages included in the JuliaFEM platform.This work aimed to implement two common model reduction methods into ModelReduction.jl:The Guyan reduction and the Craig-Bampton methods which are useful in dynamic analyses and mandatory in the commercial multi-body simulation software like AVL Excite Power Unit, which use large nite element models [49] or perform some optimization [10].As an example, Irvine has implemented these methods using Matlab [11].The ModelReduction.jlpackage is available from GitHub.
This paper introduces briey the Guyan reduction and the Craig-Bampton method and how to perform them with JuliaFEM with the help of a test model.The implemented code itself and instructions for using it are included as well.
As mentioned, FE problems widely use model reduction methods.For instance, exible multibody dynamics problems [12], problems including high-speed rotating structures [13], uid-structure interaction problems [14] and problems with localized nonlinearities, such as cracks [15], benet signicantly from model reduction.Model reduction is also useful in vibration analyses, for instance in vibration isolation modeling [16] and problems including stochastic response [17].
Structural Analyst or Researcher can apply the reduction methods to any element types, and the types and element properties can vary in the FE domain.Small and higher order elements are useful when the results change rapidly, and bigger size, loworder [18,19] elements are usable when the results are more constant [20].

Theory
Especially dynamic simulations with exible bodies require signicant computational resources.The system of equations is likely to contain a substantial number, typically of the order of millions of degrees of freedom, and require extensive computational resources to solve.To reduce the computational cost model reduction techniques are commonly used.[21,22] The model reduction methods divide into static and dynamic condensation, and dynamic condensation is a generalization of the static condensation.The following chapters present two of the most used FE model reduction techniques for both static and dynamic analyses -the Guyan reduction and the Craig-Bampton method.

Substructures and superelements
Substructuring is the process of decomposing a large FE model into smaller, componentbased models [23].It means removing elements that are unnecessary for the analysis and building larger elements -so-called superelements -out of them.These component models are called the substructures of the full system.For example, [24] views a subset of adjacent nite elements as one superelement or substructure.
Researchers use substructuring in component mode synthesis (CMS), where individual substructure problems are rst solved, and then the coupling of interfaces is built [25].CMS has many advantages in dynamic analyses especially when the assemblies are large and complex.Literature also calls substructuring and CMS as coupling problems or subsystem addition [26].One of the primary reasons for substructuring in dynamics problems is to reduce the number of degrees of freedom of the structure [27].Fewer degrees of freedom requires less computational resources than the original model.
The main steps of the substructuring process are to divide the whole structure into some substructures, to obtain reduced-order models of the components, and then to assemble a reduced-order model of the entire structure [28].Substructuring allows the evaluation of the dynamic behavior of large and complex structures.Also, local dynamic behavior can be recognized more easily by analyzing the reduced subsystems than when the entire system is analyzed [29].

Guyan reduction
Static reduction, also known as static condensation, Guyan condensation or Guyan reduction, is the most popular model reduction method presented by R.J. Guyan [30].It is a method, which ignores inertia eects of certain degrees of freedom while obtaining component modes [31].Guyan reduction is the basis for several other nite element substructuring techniques [32].
The Guyan reduction method applied in FE techniques reduces the FE model by condensing internal, dened as a slave, degrees of freedom.Specically, the technique removes the degrees of freedom located at the substructure's boundary, which is the localto-global interface.The remaining degrees of freedom dened as a master and located at the boundary, retain the stiness of the local structure but omit the inertial terms to create a denser and thus more ecient representation, at the cost of accuracy for nonstatic loading conditions.The method is only accurate for stiness reduction since Guyan reduction ignores inertial forces.[23] The static equilibrium equation is: where K is the global stiness matrix, u presents the nodal degrees of freedom and f is the nodal force vector of the static equilibrium problem.By dividing the static equilibrium equation ( 1) with regards to loaded (master) and unloaded (slave) degrees of freedom so that the forces on the unloaded degrees of freedom are zero, the static equilibrium equation becomes: where K MM , K MS , K SM , and K SS are submatrices of K and K MM is the part of K that remains after the reduction.If only u M is desired, K can be reduced as follows: where K red is the nal reduced stiness matrix.K red is obtained by writing out the set of equations as follows: Equation ( 5) can be solved for u S assuming that K SS is invertible: moreover, substituting into (4) gives Now K red can be solved as follows: where K red is the reduced stiness matrix.Structural Analyst or Researcher may choose to eliminate any component of u if the corresponding component of f is zero.The above system of linear equations is equivalent to the original equation ( 1), but it is expressed solely by the master degrees of freedom.Thus, Guyan reduction leads to a reduced system by condensing away the slave degrees of freedom.Since sparse matrix inversions require lots of computational resources, factorization methods, such as Cholesky decomposition can be applied to obtain K −1 SS and in such way reduce the calculation time.In Julia language, the Guyan reduction algorithm implementation is only a few lines of code.The code listing of the implementation is available at [33].
The Craig-Bampton method The Craig-Bampton process is a dynamic reduction technique introduced by Roy R. Craig Jr and Mervyn C. C. Bampton [34] that is widely used to assemble large-scale models (millions of degrees of freedom) that are far too computationally expensive to be modeled entirely [35].
In the Craig-Bampton process, Structural Analyst or Researcher rst separates the degrees of freedom in the original FE model into retained (master) and truncated (slave) degrees of freedom in a similar way as in the Guyan reduction.There are algorithms to help to select the master and slave degrees of freedom [36].Then, by condensing the stiness and inertial eects for the slave degrees of freedom into master degrees of freedom, the reduced model is constructed [37].
The Craig-Bampton method reduces the mass and stiness matrices of the FE model by expressing the master modes in physical coordinates and the elastic modes in modal coordinates.The method reduces the mass and stiness matrices which will contain mode shape information of the low-frequency response modes of the model.The Craig-Bampton method is especially useful in dynamic analyses that include large nite element models.[20,38] This implementation does not include damping.The equation of motion is: In the Craig-Bampton method the matrices are rst partitioned into master nodes R and slave nodes L: Equation ( 9) becomes: The division of M and K into submatrices is similar as for K in the Guyan reduction.
The degrees of freedom transform into hybrid coordinates: where I is an identity matrix, X R is a transformation matrix which relates rigid body physical displacements at the interface u R to physical displacements of the elastic degrees of freedom u L .Also, X L is a matrix of eigenvectors called normal mode shapes.It is a matrix of eigenvectors calculated from K LL and q m is a column vector of modal displacements.It is dimensionless, so all units are contained in X L .When performing the Craig-Bampton method, X L is the matrix that causes the reduction of the nal result matrices: most columns are eigenvectors to be deleted so that the matrix size reduces and the size of the nal matrices depends on the size of this matrix.Now equation ( 11) can be rewritten as To determine X R all master degrees of freedom are xed limiting consideration to a static problem (ü R = üL = 0).Equation ( 11) reduces to: The internal degrees of freedom are: where In the determination of X L the master degrees of freedom are xed.The equation of motion (9) reduces to: By assuming harmonic response (u L = X L q m e iωt ) unforced harmonic motion of the grounded structure can be expressed as: where Λ is a diagonal matrix containing the eigenvalues of (17).The eigenvectors in X L can be normalized: Since X R in ( 16) contains K −1 LL , an inverse of K LL , determining it for sparse matrices will require lots of computing resources, and it will eventually become a problem with large models.It is avoidable by determining K −1 LL as follows: Now ( 16) can be calculated as: As one can see this expression also includes an inverse, but it is an inverse of Λ which is a diagonal matrix, so it only has nonzero elements on its diagonal and therefore needs much less computing power than the computing of K −1 LL .To get the equations of motion of the system the equation ( 13) is multiplied with the transpose of the coordinate transformation matrix in (12) as follows: By simplifying, the equation of motion (9) becomes Above is the nal form of the dynamic equation of motion for the Craig-Bampton method when the generalized mass matrix is normalized, damping is ignored, and only boundary (master) forces are in consideration, which is true for most practical problems [38].For the JuliaFEM implementation equation ( 24) is expressed as: where The code listing of the Craig-Bampton method JuliaFEM implementation is available at [33].

Test model
The JuliaFEM model reduction algorithms are tested on an example model to verify the codes work correctly.
The example model is a 1-dimensional rod with four elements and ve nodes.The rod is xed at node 1, and it also has four roller supports at nodes 2 -5.The rollers are not necessary since the model is a rod, but they included in the model since in some commercial FEM programs the rod dividing nodes are interpreted as joints, so that horizontal support is needed when performing the dynamic analysis.Because of these supports, node 1 has 0 degrees of freedom and nodes 2 -5 have 1 degree of freedom.There is a horizontal driving force at node 5.The model is presented in Figure 1 where the length of one element is L = 0.25 and the driving force at node 5 is F = 1 N.
It is evident that the model does not need as many elements.Since the model is a rod, only one element would be enough to give the correct displacement at node 5.The model reduction will reduce the mesh so that only one element is left.
Before the model reduction, authors calculated the standard static and modal analyses.Then the model reduction methods were performed with the implemented functions, and the authors performed analyses with the reduced matrices.Finally, authors compared the results of both analyses.

Static analysis
For the example model, equation ( 1) is the following without considering boundary conditions: where the u-vector is the solution of the equation ( 1).Only the displacements u 1 and u 5 are globally meaningful.The static condensation will remove the undesired degrees of freedom and give the same result with much smaller matrices.

Modal analysis
The global stiness and mass matrices for the example model, when the boundary conditions are taken to account, are the following: The equation of motion (9) will give the eigenvalue problem: yielding the eigenmodes x of the model and the eigenvalues whose square roots are the angular eigenfrequencies ω.For the example model the eigenvalues calculated from (31) are the following: The calculation of natural frequencies is as follows: where f n are natural frequencies.Equation (33) gives the following frequencies for the example model: The dynamic reduction will condense the stiness and mass matrices of the model and give fewer eigenmodes, but the new modes are among the original low-frequency response modes.
Reduction methods applied to the example model Now the Guyan reduction and the Craig-Bampton method will be applied to the example model.The substructuring will remove nodes 2 -4 and there will be only one element left the superelement.

Guyan reduction by hand
For the example model, the submatrices in (2) will be the following since only u 1 and u 5 are desired degrees of freedom.Figure 3 shows how the model's stiness matrix ( 26) is divided into submatrices by the desired degrees of freedom.
Figure 3: Dividing K into submatrices.Now equation ( 8) becomes the following: which gives The equation ( 3) with K red being given by can now be solved for u and f to obtain u 1 = 0 and u 5 = 4.The result is in agreement with the solution vector in (28) for the chosen master unknowns.where K, M, r and l and n are original stiness matrix, original mass matrix, master degrees of freedom, slave degrees of freedom and the number of the internal modes to keep respectively.Users may choose r, and l the way they wish and n so that n ≤ length of l, remembering that the size of these variables will aect the size of the result matrices so that a small n gives small matrices.When the subpackage ModelReduction.jl of the JuliaFEM platform is used, the Craig-Bampton method is applied by simply calling the Craig-Bampton function craig_bampton(), which gives two matrices as an output.Mred is the reduced mass matrix and Kred is the reduced stiness matrix of the model.The sizes of the reduced matrices are (r+n)x(r+n).The number of modes that is to be computed from the reduced matrices is r+n.
Table 1 shows the natural frequencies of the example model computed with the reduced matrices with dierent values of n compared to the frequencies computed with the original stiness and mass matrices.1 shows that the dierence between the frequencies calculated with the reduced matrices and the original matrices increases when n decreases.The quality of the eigenvalue approximation will decrease as the mode number increases.

Large examples
In the following the Craig-Bampton method in ModelReduction.jl is performed on large scale models.The rst example is is a 3D-model of a bracket that is attached to two adapter plates via tie contacts and xed from the plates.The original model has almost 300 000 degrees of freedom, but after performing the reduction, only about 200 degrees of freedom remains.Figure 4 presents the model and the locations of the master degrees of freedom.
Natural frequencies have been calculated before with JuliaFEM to this example [3].Table 2 shows the ve lowest frequencies calculated with the reduced matrices compared to the original frequencies and the dierence between them in percent when the number of master degrees of freedom is 192 and using ten of the generalized eigenvalues.The dimensions of the stiness and mass matrices reduce from 293310 to 202. Figure 5 presents the rst eigenmode of the bracket.The source code in Julia syntax for the model reduction and natural frequency analysis as shown in listing 1.This paper presents the Guyan and Craig-Bampton methods implementations in the Ju-liaFEM open source platform for nite element modeling development.Also, it demonstrates how compelling the Julia programming language is for scientic simulations as well as how robust the JuliaFEM platform is for nite element method development.Just a few lines of code implements a new missing feature, and the correct use of Julia will guarantee a speed close to that allowed by the C-language.
The condensed stiness and mass matrices give the same results as the original matrices which prove that the implemented algorithms work correctly.The inaccuracy of the dynamic condensation increases when the number of retained internal modes decreases, and it also depends on selected master degrees of freedom.The lowest frequencies calculated with the condensed matrices are the most accurate.

Figure 1 :
Figure 1: The original test model.

Figure 2
presents the new structure.The new variables of the model are the same except the length L of the element since it now refers to the length L = 1.0 of the whole rod.The following paragraphs present detailed steps of the process leading to this superelement.
Guyan reduction applied to the example modelInstructions regarding the Guyan reduction function, included in the ModelReduction.jlsubpackage of the JuliaFEM platform, are presented next.The following example of using the function is Julia syntax.The rst step is to install ModelReduction.jlpackage.Then the Guyan reduction function usage is as follows: julia> Pkg.add("ModelReduction") m, and s are original stiness matrix, master nodes, and slave nodes respectively.When Structural Analyst or Researcher uses the subpackage ModelReduction.jl of the JuliaFEM platform, the Guyan reduction is applied by simply calling the Guyan reduction function guyan_reduction(), which gives one reduced matrix as a result.Kred is the reduced stiness matrix of the model.The Craig-Bampton method applied to the example model Next example demonstrates how to use the Craig-Bampton function included in the Mod-elReduction.jlsubpackage of the JuliaFEM platform.The following example is Julia syntax and demonstrates the use of reduction function to the example model.First, the ModelReduction.jlpackage must be installed (if it is not yet the case) and the following variables need to be dened: julia> Pkg.add("ModelReduction") Kred = ModelReduction.craig_bampton(K,M, r, l, n) ([2.75 -1.20711; -1.20711 1.0], [0.25 0.0; 0.0 0.292893]),

Figure 4 :
Figure 4: The bracket model with nearly 300 000 degrees of freedom and the locations from where the master degrees of freedom were chosen [3].

Figure 6 :
Figure 6: The industrial sized model with 12.6 million degrees of freedom and the locations from where the master degrees of freedom were chosen [2].

Figure 7 :
Figure 7: The rst eigenmode of the industrial sized model .

Table 1 :
The natural frequencies computed with the original K and M compared to the frequencies computed with Kred and Mred with dierent sizes of n.

Table 2 :
Natural frequencies of the bracket compared to the natural frequencies of the reduced model.Mode Original f [Hz] Reduced f [Hz] Dierence [%]

Table 3 :
The dierence between frequencies of the industrial size model compared to the natural frequencies of the reduced model in percent.