Study of the model order reduction of an entire aircraft configuration

Secci´o de Terrassa Centre: ETSEIAT Lecture: Departament de Resist`encia de Materialas Work Title: Final Master Thesis Study of the Model Order Reduction of an Entire Aircraft Configuration Grupo: M2B-P2016/17 Delivery Date: 20th January, 2017 Students: Fabi´an Lajas Contreras Director: Prof. Joaqu´ın Hern´andez Ortega CONTENTS Contents Contents 2 List of Figures 3 List of Tables 6 1 Introduction 9 1.1 AimoftheStudy..................................... 9 1.2 State-Of-The-Art..................................... 9 1.2.1 Origins of Substructuring Methods . . . . . . . . . . . . . . . . . . . . . . . 9 1.2.2 Reduced Order Models based on substructuring . . . . . . . . . . . . . . . . 10 1.3 Scope ........................................... 10 1.4 Requirements....................................... 10 1.5 Employedtools...................................... 11 1.6 Outlineofthestudy ................................... 11 2 Theoretical Basis 11 2.1 ModalAnalysis...................................... 11 2.2 BasicFEMconcepts ................................... 13 2.2.1 The 2-Noded Euler Bernouilli Beam Element . . . . . . . . . . . . . . . . . 13 2.2.2 Reissner-Mindlin Flat Shell Element . . . . . . . . . . . . . . . . . . . . . . 15 2.3 Substructuring Techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3.1 Craig-Bampton Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3 Application for Beam Geometry 19 3.1 BasicCodeStructure................................... 19 3.2 BeamProblemDefinition ................................ 21 3.3 Results Extracted for a Beam Problem . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.4 Discussionofresults ................................... 26 4 Application for a Shell Geometry 26 4.1 Editing Kratos ..................................... 28 4.2 ModifiedCodeStructure................................. 31 4.3 Reissner-Mindlin Flat Shell Element Problem Definition . . . . . . . . . . . . . . . 33 4.4 Results extracted for a Reissner-Mindlin Flat Shell Element . . . . . . . . . . . . 35 4.5 Discussionofresults ................................... 39 5 Application for a Complete Aircraft Configuration 39 ETSEIAT. Fabi´an Lajas Contreras 2 5.1 ChosingtheGeometry.................................. 39 5.1.1 KindsofModels ................................. 40 5.1.2 Chosenmodel................................... 41 5.2 FinalCodeStructure................................... 42 5.3 Complete Airplane Problem Definition . . . . . . . . . . . . . . . . . . . . . . . . . 43 5.3.1 DesignProcess .................................. 43 5.3.2 DesignProblems ................................. 47 5.3.3 Final Geometry Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 5.4 Results .......................................... 48 5.5 Discussion......................................... 54 6 Concluding remarks 54 7 Future Lines of Research 55 References 56 A Costs 57 A.1 WorkingsHours...................................... 57 A.2 SoftwareLicenses..................................... 58 A.3 Hardware ......................................... 58 A.4 TotalBudget ....................................... 58 B Environmental Impact 59 List of Figures 1 Example of a substructured system; the complete structure has been broken down into six different components. Image courtesy of [1]. . . . . . . . . . . . . . . . . 10 2 Two-noded Euler-Bernoulli beam element, image couurtesy of [2]. . . . . . . . . . 14 3 Mesh discretization, it can be distinguished two kinds of shells, the quadrilateral mesh at left, and the triangular mesh at right. For this study, the triangular element will be used. Image courtesy of [2]. . . . . . . . . . . . . . . . . . . . . . . . . . . 15 4 Partitioning of a Structure into Two Substructures. In the figure Γ refers to the boundary between two substructures. Image courtesy of [3]. . . . . . . . . . . . . 16 5 This is the first iteration of the code. . . . . . . . . . . . . . . . . . . . . . . . . . 21 6 Fokker D.VII blueprints, it can be aprreciated the truss structure. . . . . . . . . . 22 7 Total model, the structure nodes appear in black. . . . . . . . . . . . . . . . . . . 22 8 The restricted nodes for this modal analysis will be located at the bottom of the structure.......................................... 23 ETSEIAT. Fabi´an Lajas Contreras 3 LIST OF FIGURES 9 Complete model. Each substructure has been marked with a different colour: Green for the first substructure, Yellow for the second substructure, Blue for the third substructure........................................ 23 10 Modal Shapes Obtained for each substructure. It is important to comment that in these analysis, the boundary DOF’s are considered as imposed DOFs. . . . . . . . 24 11 First Modal Shape, a) refers to the classic algorithm: CA, b) refers to the substructuredalgorithm:CB. .................................. 24 12 Second Modal Shape, a) refers to the classic algorithm: CA, b) refers to the substructuredalgorithm:CB. ............................... 25 13 Third Modal Shape, a) refers to the classic algorithm: CA, b) refers to the substructuredalgorithm:CB. ............................... 25 14 Fourth Modal Shape, a) refers to the classic algorithm: CA, b) refers to the substructuredalgorithm:CB. ............................... 26 15 In this diagram it is shown the structure of a Kratos simulation. . . . . . . . . . 27 16 Basic input file mdpa, it can be appreciated the nodes tag, the connectivity tag and the conditions tag. Image courtesy of [4]. . . . . . . . . . . . . . . . . . . . . 28 17 Applications folder inside the main Kratos directory. In the upper part of the window, it can be seen the installation directory, whereas the arrow, also in green is pointing towards the SolidMechanicsApplications folder. . . . . . . . . . . . . . 29 18 Location of the file add custom utilities to python.cpp inside Kratos directory. . 30 19 Location of the file print matrix.h inside Kratos directory. . . . . . . . . . . . . 30 20 Application directory. The name of the application launcher is marked with green MainKratos.py. ..................................... 31 21 Substructural ALgorithm Diagram, including Kratos module. . . . . . . . . . . . 32 22 P51 mustang, chosen airplane for the second problem. . . . . . . . . . . . . . . . . 33 23 Image Base imported into Libre CAD a). Points Extracted using the same software, b) later this points are converted into .svg format, which can be read using Matlab TM. The final geometry is depicted in c). .............. 34 24 Restricted nodes chosen for this problem are marked in red —the ones with lowest Zcoordinate........................................ 34 25 Substructured model chosen for this problem. Each color is identified with a differentsubstructure. ................................... 35 26 Modal Shapes Obtained for the first and second substructures. Recall that boundary DOF’s are considered as imposed DOFs. . . . . . . . . . . . . . . . . . . . . . 35 27 Modal Shapes Obtained for the third and fourth substructures. . . . . . . . . . . 36 28 Modal Shapes Obtained for the fifth and sixth substructures. . . . . . . . . . . . 36 29 First Modal Shape, a) refers to the classic algorithm: CA, b) refers to the substructuredalgorithm:CB. .................................. 37 30 Second Modal Shape, a) refers to the classic algorithm: CA, b) refers to the substructuredalgorithm:CB. ............................... 38 ETSEIAT. Fabi´an Lajas Contreras 4 31 Third Modal Shape, a) refers to the classic algorithm: CA, b) refers to the substructuredalgorithm:CB. ............................... 38 32 Fourth Modal Shape, a) refers to the classic algorithm: CA, b) refers to the substructuredalgorithm:CB. ............................... 39 33 First kind of airplanes, classic models A), side view B), profile view C) top view . 40 34 Second kind of airplanes, militar models A), side view B), profile view C) top view 40 35 Third kind of airplanes, private aviation A), side view B), profile view C) top view 40 36 Cessna Citation. More information of this product can be found in [5]. . . . . . . 41 37 Finalcodestructure................................... 42 38 Airplane scheme. Five main parts can be distinguished: wings in Green, engine structure and fairing in Blue, tail in Pink, Stabilizers in Grey. ........... 43 39 Wing structure: on the RIGHT, scheme structure extracted from the Airbus A-220 guide: on the LEFT, designed wing structure. Notice that no control surface is included. ......................................... 44 40 Scheme of the fuselage approximation into a cylinder: a) Scheme of the stringer creation used in the geometry design. b) Stringers Cross Sections in BLACK, Guide Curves in RED, Frames in BLUE. . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 41 Stringer sections designed: a) Original section, based upon a common stringer profile extracted from a standard airliner. b) Chosen Approximation (T section with pointing upper side). c) Extruded shape applied in the model, with the union between the stringer and the frames. d) Cross section T. . . . . . . . . . . . . . . 45 42 a) Profile of the fuselage, b) Front view of the fuselage, c) Isometric view. . . . . . 45 43 Horizontal Tail Structure: on the LEFT scheme structure extracted from the Airbus A-220 guide, on the RIGHT, designed tail structure (control surfaces have been neglected). ........................................ 46 44 Vertical Stabilizer: on the LEFT, scheme structure extracted from the Airbus A220 guide, on the RIGHT, designed vertical stabilizer structure, control surfaces such as the rudder have been neglected. . . . . . . . . . . . . . . . . . . . . . . . . 46 45 Fairing Structure and Pylon: on the LEFT, schematic structure of the pylons extracted from the Airbus A-220. On the RIGHT, representation of the model created with Solid Works TM............................. 46 46 Union between the wing and the fuselage. Mesh compatibility were achieved after several iterations (geometry modifications). . . . . . . . . . . . . . . . . . . . . . 47 47 Main views of the final model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 48 Dirichlet boundary conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 49 Structurepartition. ................................... 48 50 Modal Shapes for the first, and the second substructures. . . . . . . . . . . . . . 49 51 Modal Shapes for the third and fourth substructures. . . . . . . . . . . . . . . . . 50 52 Modal Shapes for the fifth and sixth substructures. . . . . . . . . . . . . . . . . . 51 53 FirstModalShape.................................... 52 54 SecondModalShape .................................. 53 ETSEIAT. Fabi´an Lajas Contreras 5 LIST OF TABLES 55 ThirdModalShape ................................... 53 56 FourthModalShape .................................. 54 57 Mesh refinement: a) 17005 nodes with 36651 elements, b) 35599 nodes with 75739 elements, c) 128794 nodes with 268336 elements. . . . . . . . . . . . . . . . . . . . 55 List of Tables 1 Materials Table Chosen for this first problem. . . . . . . . . . . . . . . . . . . . . . 22 2 In this table, it can be seen the modal analysis performed at each substructure. . . 23 3 Natural Frequecies, CA refers to classical algorithm , whereas CB refers to CraigBampton algorithm algorithm. ........................... 24 4 Reissner-Mindlin Flat Shell Element properties . . . . . . . . . . . . . . . . . . . . 34 5 Modal frequencies of each substructure. . . . . . . . . . . . . . . . . . . . . . . . . 37 6 Natural Frecuecies for the shell problem, CA refers to Classic Algorithm, whereas CB refers to Craig-Bampton algorithm. . . . . . . . . . . . . . . . . . . . . . . . . 37 7 Cessna Citations XLS Specifications, data courtesy of [6] . . . . . . . . . . . . . . 41 8 Materialproperties. ................................... 48 9 Natural frequencies for each substructure . . . . . . . . . . . . . . . . . . . . . . . 52 10 Natural frequencies. CA refers to Classic Algorithm, whereas CB refers to the Craig-Bamptonalgorithm................................. 52 11 Workinghourscosts ................................... 57 12 Softwarelicensescosts .................................. 58 13 Hardwarecosts ...................................... 58 14 Workinghourscosts ................................... 58 ETSEIAT. Fabi´an Lajas Contreras 6 Acknowledgements I would like to thank to everyone who has helped me in the planification, development and achievement of this final master project, namely, to •Joaqu´ın A. Hern´andez Ortega •Riccardo Rossi •My parents and siblings ETSEIAT. Fabi´an Lajas Contreras 7 LIST OF TABLES Nomenclature ¨ dAcceleration. A B Array of constants. DDumping Matrix. dDisplacement. λiEigenvalue. FForce Vector. IIdentity Matrix. KStiffness Matrix, submodel Stiffness Matrix. KBB Boundary Stiffness Matrix. KBI Interface Stiffness Matrix. KII Internal Stiffness Matrix. XIReduced Modal Shape. MMass Matrix, submodel Mass Matrix. RReal Numbers. MBB Boundary Mass Matrix. MBI Interface Mass Matrix. MII Internal Mass Matrix. φModal Vector. ΩSpectral Matrix. ωNatural Frequency. pExternal Forces. qAmplitude. RCBT M Craig-Bampton Transformation Matrix. ˜ KReduced Stiffness Matrix. ˜ MReduced Mass Matrix. ρMaterial density. ˜ KBB Reduced Inner Stiffness Matrix. ˜ KTTotal reduced Stiffness Matrix. ˜ MBReduced Inner Mass Matrix. ˜ MTTotal reduced Mass Matrix. RCBT M Reduced Craig-Bampton Transformation Matrix. uBBoundary Results uCCompressed Results uTotal Results ˙ dVelocity. xTotal degrees of Freedom. xBSubsystem internal degrees of freedom. xIBoundary degrees of freedom. 0Array or matrix composed by zeros. GTorsional Inertia Product. IxInertia Product in the x axis. IyInertia Product in the y axis. IzInertia Product in the z axis. lLength of a slender beam. mNumber of reduced DOFs. nNumber of elements present in an enumeration. nBNumber of boundary DOFs. tShell Thickness. A Cross sectional area for a slender beam. E Young Modulus. t Time. ETSEIAT. Fabi´an Lajas Contreras 8 1 Introduction 1.1 Aim of the Study The primary aim of this study is to understand and implement in a computer program the basic theory underlying model order reduction of Finite Element (FE) structural models. More concretely, attention is focused on modal analysis via substructuring techniques, and the interest lies in comparing the performance of such substructuring methods with the standard modal analysis of linear elastodynamic FE models. The method is applied to 3 aeronautical structures, ranging in complexity from a simple truss structure to a complete airframe of a private jet aircraft. For the truss structure, all required operations are carried out with a Matlab code developed by the author. For the other two structures, the finite element information ( mass and stiffness matrices ) is retrieved from an open source FE software called KRATOS, and then processed by a Matlab script in order to determine the vibrational modes and the natural frequencies. It should be highlighted that the construction of the geometric models of the three tested airframes have been also developed, from scratch, by the author. 1.2 State-Of-The-Art 1.2.1 Origins of Substructuring Methods Substructuring methods were invented in the early 1960s by aircraft engineers to carry out a firstlevel breakdown of complex systems such as a complete aircraft shown in Figure 1. The original factors that motivated the development of this kind of techniques were, according to [1]: •To Facilitate division of labor: Substructures with different functions are done by different teams. This fact would allow different groups to work separately with different parts of the project, this fact would allow to perform geometry modifications, nevertheless it is important to remind that the interface nodes that each substructure contains must remain the same. •Take advantage of repetition: It is usual for this kind of projects to have symmetric parts in the structure, for instance, it is common to deal with aircrafts which contain symmetry planes that split by the half the entire structure, recognizing repetitions, saves model creation and, of course, time. •Overcome computer limitations: This kind of techniques allow to analyze entire huge models part by part, saving computational time and surpassing the memory limitations that a usual machine may have. ETSEIAT. Fabi´an Lajas Contreras 9 2 THEORETICAL BASIS 4. Linear assumed transverse shear strain field is assumed as in other Reissner-Mindlin element formulation: the TLQL. More information about this shell model can be found in [2]. 2.3 Substructuring Techniques Substructuring methods are frequently used in dynamic analysis; the main reasons that motivated the development of these methods are: •The interest in lower frequency eigensolutions, which is common in structural analysis, and it may prove to be advantageous to reduce the starting problem to a smaller dimension. •In large projects, such as aircrafts, separate parts of the complete system can be analized by each project team, later these parts will be used to reconstruct the model of the whole system. Therefore, the interconnecting surfaces must be defined carefully in order to ensure the compatibility of the different parts of the model, as it can be seen in Figure 4. Figure 4: Partitioning of a Structure into Two Substructures. In the figure Γ refers to the boundary between two substructures. Image courtesy of [3]. 2.3.1 Craig-Bampton Method The Craig-Bampton method, (Craig and Bampton 1968), also commonly referred to as the component modes method or the modal synthesis is mainly used as a dynamic substructuring method. It consist in reducing the complexity —via standard modal analysis— of each part of the model. The first step in the formulation of this method is the separation into inner nodes and boundary nodes: x=  xI xB  ,K=  KII KIB KBI KBB  ,M=  MII MIB MBI MBB  and F=  0 pB  (20) Considering the modal analysis equation, Eq. (21), the total system can be rewritten by using the previously defined equations; the final result is Eq. (22). Kx−ω2Mx=F(21)   KII KIB KBI KBB    xI xB  −ω2  MII MIB MBI MBB    xI xB  =  0 pB  (22) ETSEIAT. Fabi´an Lajas Contreras 16 Considering that if there were no inertia forces, the internal degrees of freedom could be computed by static condensation, it is possible to define the xIas: xI=−K−1 II KIBxB(23) Boundary DOFs xBwill be considered as static DOFs, and an eigenvaule problem will be solved by considering the internal DOFs as the free ones: KII XI=ω2 IMII XI(24) Of course, according to the modal problem, the fixed interface modes are stored in the columns of matrix XIso that: XIKII XI=diag(ω2 I1... ω2 I nI ) = Ω2 I(25) XIKII XI=I(26) It is possible to build the total DOFs x by using the following relation: x=RCBT M   η xB  (27) where RCBT M is what it is called the Craig Bampton Transformation Matrix: RCBT M =  XI−K−1 II KIB 0I (28) and ηrefers to the intensity parameters of the substructure’s internal vibration modes. Once all these elements have been defined, the initial model of each substructure will be reduced. Only a certain number of internal modes Iwill be chosen: XI=hxI(1) ... xI(m)i(29) On the other hand, in this reduction process, it is necessary to maintain the boundary DOFs in order to assure deformation compatibility at the interfaces; the final form of CBTM will be: RCBT M =  XI−K−1 II KIB 0I (30) The dimension of this matrix will be n×(m+nB). The reduced stiffness and mass matrix will be extracted by performing the following operations: ˜ K=RT CBT M KRCBT M and ˜ M=RT CBT M MRCBT M (31) For each submodel: ETSEIAT. Fabi´an Lajas Contreras 17 2 THEORETICAL BASIS ˜ M=  ˜ Ω2 I0 0˜ KBB  and ˜ K=  I˜ MB ˜ MB˜ KBB  (32) If a problem with N substructures is considered, every substructure will have a transformed stiffness matrix and a transformed mass matrix. For the i-th substructure, it is defined as ˜ Kiand ˜ Mi. The total number of boundaries present in the model can be defined as M. Then the final assembled matrix for the entire system can be defined as: ˜ KTand ˜ MT: ˜ KT=  ˜ Ω2 I{1...m}0 0˜ KBB {1...nB}  (33) ˜ MT=  I{1...m}˜ MBB {1...m,1...nB} ˜ MBB {1...nB,1...m}˜ MBB {1...nB}  (34) These matrix will have the following expanded form: ˜ KT=               ˜ Ω2 I 1 . . . 0 0 . . . 0 . . ..... . .. . ..... . . 0. . . ˜ Ω2 I m 0. . . 0 0. . . 0˜ KBB 1+˜ KBB j. . . 0 . . ..... . .. . ..... . . 0. . . 0 0 . . . ˜ KBB 1+˜ KBB j               (35) ˜ MT=               I1. . . 0˜ MB 1 . . . 0 . . .Ij . . .˜ MB j ...˜ MBnB 0. . . Im0. . . ˜ MBm, nB ˜ MB 1 ˜ MB j 0˜ MBB 1+˜ MBB j. . . 0 . . ..... . .. . ..... . . 0˜ MBnB˜ MBnB, m 0. . . ˜ MBB 1+˜ MBB j               (36) Now that the entire system is assembled, the required analysis is applied. Then the modal analysis is employed once the entire system is assembled. After that, it will be necessary to decompress the obtained modal vectors . This process can be done by using the equation Figure 37. In this equation uCrefers to the compressed results, uBdenotes the boundary results, and finally uare the final decompressed results of the analysis. hui=RCBT M   uC uB  (37) ETSEIAT. Fabi´an Lajas Contreras 18 3 Application for Beam Geometry In this part, the main code structure required to solve a modal analysis using a substructured technique will be introduced. It will be used Matlab TM in order to create the main algorithm. On the other hand, GiD TM will be used for the meshing process. In the following section the code structure will be deeply studied. Later, a beam problem will be solved and the results obtained will be compared and commented. 3.1 Basic Code Structure The program will compute the modal analysis of an aeronautic structure: the fuselage, the wings, the stabilizers, the landing gears, or even entire body configurations. As it was said in Section ?? , the results obtained by two methods will be compared. It is for this reason that the created code will contain the following algorithms: •the Classic modal analysis algorithm, which consists in the analysis of the entire geometry without using any kind of substructuring method. This process may prove slow, specially for large geometries. Hereafter, we shall refer to this process as classical algorithm . •The Craig-Bampton modal analysis algorithm. This algorithm will analyze the entire geometry using the substructuring method. Hereafter, it will be referred as Craig-Bampton algorithm . The code structure comprises several functions and subfunctions created using the vectorization technique. According to [14], vectorization consists in the transformation of a code in order to use nonscalar objects. Using this technique, speedup factors of almost 10 can be reached. In Figure 5, the first code structure is presented, which is composed by the following kinds of objects: 1. ORANGE TAGS are used as starting and ending execution points, and do not refer to functions. For this kind of tag the following parts are defined: (a) START: This part is mainly used as the starting point of the program. (b) END: As it was commented in the last item, this tag works as a state of the code. In this case it works as the ending point of the program. 2. BLACK TAGS refer to decision points. These decision points are triggered by a boolean object that must be correctly filled before starting the program execution. (a) EDITMODE: Refers to a conditional state. Depending on the value of the boolean object, it will take one or another path. (b) SUBSTRUC:EDITMODE will behave in the same way as the previous tag. 3. RED TAGS are used to mark functions that perform an auxiliary functionality. (a) VALIDATOR: Validator refers to an external functionality that assures that all the inputs present in the program are correctly introduced. 4. BLUE TAGS refer to the main functionalities run by the program. (a) INPUT SETTINGS: In this part, the main settings of the program are defined. ETSEIAT. Fabi´an Lajas Contreras 19 3 APPLICATION FOR BEAM GEOMETRY (b) INPUT GEOMETRY RECIPE: In this part, the geometry to analyze is defined. (c) CREATE GEOMETRY : This function is in charge of creating the geometry using the information introduced in the previous state. (d) READMSH : In this part, the mesh of the problem and the main conditions (materials and fixed DOFs), are extracted. (e) MATRIX ASSEMBLY : This part is only activated if the SUBSTRUC switch is turnedoff, it is in charge of computing the elemental stiffness and mass matrix, perform the assembly and applying the conditions defined in READMSH. (f) SUBSTRUC ALGORITHM : This part is activated only if the SUBSTRUC switch is turned-on. This part will divide the mesh components into different substructures; it also will dump all the information into the hard drive in order to save memory. (g) MATRIX ASSEMBLY SUB: As in the previous tag, SUBSTRUC ALGORITHM, this part can only be accessed by switching-on SUBSTRUC ; it will read each one of the substructures, perform the Craig-Bampton transformation, and finally, compute the final assembly matrix. (h) EIGENVAL ALGORITHM : This part will solve an eivenvalue problem. (i) MODAL SHAPE ALGORITHM : This part will compute the modal shapes extracted by the previous tag. EIGENVAL ALGORITHM. 5. GREEN TAGS This kind of tag refers to the external software that it is used to perform special functionalities, such as the meshing algorithm, that is run by GiD TM and the matrix computing and printing system, which is done by Kratos. (a) GID MESHER: External functionality can only be accessed by switching-off EDITMODE; it calls the software GiD TM in order to load the created geometry and perform the meshing process. (b) GID VIEWER: External functionality that can only be accessed by switching-on EDITMODE; it calls the software GiD TM in order to visualize the created geometry. ETSEIAT. Fabi´an Lajas Contreras 20 Figure 5: This is the first iteration of the code. 3.2 Beam Problem Definition In this section, a modal analysis using the Craig-Bampton algorithm is performed. A 3D beam structure will be tested using the models defined in Section 2.2.1 . The beam element is a reasonable solution for truss-like geometries such as the helicopter fuselage or classic airplanes such as those shown in Figure 6. This part will serve as a test for the Craig-Bampton algorithm . ETSEIAT. Fabi´an Lajas Contreras 21 3 APPLICATION FOR BEAM GEOMETRY Figure 6: Fokker D.VII blueprints, it can be aprreciated the truss structure. The geometry consists in a airplane-like structure composed by beams, as it can be seen in Figure 7. The model is formed by circular cross-section beams made of the same material. The physical properties chosen for this element are displayed in Table 1. Figure 7: Total model, the structure nodes appear in black. Table 1: Materials Table Chosen for this first problem. E[MPa]ν G [MPa]Rb[m]Ab[m2]Ix[m4]Iy[m4]Iz[m4]ρ[Mkg/m3] 206900 0.29 80193.7984 0.025 0.196e-2 6.136e-007 3.068e-007 3.068e-007 7.85e-003 ETSEIAT. Fabi´an Lajas Contreras 22 Figure 8: The restricted nodes for this modal analysis will be located at the bottom of the structure. Just before solving the problem all the substructures are defined. In Figure 9 the chosen substructures for this example are presented. Figure 9: Complete model. Each substructure has been marked with a different colour: Green for the first substructure, Yellow for the second substructure, Blue for the third substructure. 3.3 Results Extracted for a Beam Problem The number of inner DOFs have been set to 6. In Table 2, the natural frequencies for each one of the substructure are presented. The three lowest modal shapes for each substructure are shown in Figure 10. On the other hand, the tenth lowest natural frequencies for the total analysis can be found in Table 3. Finally the 4 lowest modal shapes are presented in Figure 11, Figure 12, Figure 13 and Figure 14. Table 2: In this table, it can be seen the modal analysis performed at each substructure. Freq [Hz] ω1ω2ω3ω4ω5ω6 sub 1 2.2712 9.8239 14.2458 14.5011 15.5774 16.6988 sub 2 22.1827 52.4582 54.3966 81.5935 113.3861 124.2681 sub 3 2.4120 10.0043 14.4276 15.3735 16.2497 17.0245 ETSEIAT. Fabi´an Lajas Contreras 23 3 APPLICATION FOR BEAM GEOMETRY Figure 10: Modal Shapes Obtained for each substructure. It is important to comment that in these analysis, the boundary DOF’s are considered as imposed DOFs. Table 3: Natural Frequecies, CA refers to classical algorithm , whereas CB refers to CraigBampton algorithm algorithm. Freq [Hz]ω1ω2ω3ω4ω5ω6ω7ω8ω9ω10 CA 1.651 1.684 8.159 8.563 8.741 13.953 15.548 15.677 24.091 24.104 CB 1.647 1.676 8.530 8.550 11.791 15.532 15.647 18.330 23.958 24.089 Figure 11: First Modal Shape, a) refers to the classic algorithm: CA, b) refers to the substructured algorithm: CB. ETSEIAT. Fabi´an Lajas Contreras 24 Figure 12: Second Modal Shape, a) refers to the classic algorithm: CA, b) refers to the substructured algorithm: CB. Figure 13: Third Modal Shape, a) refers to the classic algorithm: CA, b) refers to the substructured algorithm: CB. ETSEIAT. Fabi´an Lajas Contreras 25 4 APPLICATION FOR A SHELL GEOMETRY Figure 21: Substructural ALgorithm Diagram, including Kratos module. ETSEIAT. Fabi´an Lajas Contreras 32 4.3 Reissner-Mindlin Flat Shell Element Problem Definition For this problem the same analysis that was done in Section 3.2 will be performed. The chosen geometry corresponds to an approximation of the actual geometry of a Mustang P-51, displayed in Figure 22. Several parts have been removed, such as the propeller, the landing gear, the cockpit and the control surfaces. This decision was make in order to save editing time and computational effort. Figure 22: P51 mustang, chosen airplane for the second problem. The geometric model has been created from scratch by using a GNU software called Libre CAD . Several representative points of the surface of the plane were extracted from the picture shown in Figure 22 by using such a software. The coordinates of such points were then imported into Matlab TM by using a homemade script, and the final shape was drawn from such points using simple geometric transformations. The final result can be seen in Figure 23. For this problem, Reissner-Mindlin Flat Shell Element is also used. Based on the recommendation given in [15], [16] and [17], an Aluminium 2014 has been chosen; its properties are displayed in Table 4. It should be mentioned that the wings have been reduced into plates without considering nor the wings thickness neither the airfoils. Internal structure such as frames and stringer have been also washed out. To compensate for this lack of internal structure, a relatively high thickness has been set for all shell elements (t= 23[mm]). The restricted nodes will be defined similarly as in the previous problem, see Figure 24. Finally, the employed substructured model is depicted in Figure 25. ETSEIAT. Fabi´an Lajas Contreras 33 4 APPLICATION FOR A SHELL GEOMETRY Figure 23: Image Base imported into Libre CAD a). Points Extracted using the same software, b) later this points are converted into .svg format, which can be read using Matlab TM. The final geometry is depicted in c). Table 4: Reissner-Mindlin Flat Shell Element properties E[MPa]ν t [mm]ρ[kg/m3] 185 0.29 23.62 2795.6704 Figure 24: Restricted nodes chosen for this problem are marked in red —the ones with lowest Z coordinate. ETSEIAT. Fabi´an Lajas Contreras 34 Figure 25: Substructured model chosen for this problem. Each color is identified with a different substructure. 4.4 Results extracted for a Reissner-Mindlin Flat Shell Element For this problem the inner DOFs have been reduced to 6, in Table 5, the natural frequencies for each defined substructure are presented. Also the 3 lowest modal shapes for each substructure are presented in Figure 26, Figure 27 and Figure 28. The tenth lowest natural frequencies for the total analysis are presented in Table 6, while the 4 lowest modal shapes are presented in Figure 29, Figure 30, Figure 31 and Figure 32. Figure 26: Modal Shapes Obtained for the first and second substructures. Recall that boundary DOF’s are considered as imposed DOFs. ETSEIAT. Fabi´an Lajas Contreras 35 4 APPLICATION FOR A SHELL GEOMETRY Figure 27: Modal Shapes Obtained for the third and fourth substructures. Figure 28: Modal Shapes Obtained for the fifth and sixth substructures. ETSEIAT. Fabi´an Lajas Contreras 36 Table 5: Modal frequencies of each substructure. Freq [Hz] ω1ω2ω3ω4ω5ω6 sub 1 1.048260 1.713267 1.795929 2.051887 2.119760 2.142968 sub 2 0.848251 1.967279 2.096361 2.837432 3.129295 3.348210 sub 3 0.692263 1.840554 2.112059 2.724539 2.945360 3.481899 sub 4 0.240465 0.627079 1.019696 1.034580 1.311245 1.427693 sub 5 0.237343 0.628036 1.016853 1.027050 1.300488 1.430072 sub 6 0.692263 1.840554 2.112059 2.724539 2.945360 3.481899 Table 6: Natural Frecuecies for the shell problem, CA refers to Classic Algorithm, whereas CB refers to Craig-Bampton algorithm. Freq [Hz]ω1ω2ω3ω4ω5ω6ω7ω8ω9ω10 CA 0.2306 0.2347 0.3720 0.5363 0.6200 0.6335 0.6873 0.7247 0.7341 0.8455 CB 0.2306 0.2347 0.3732 0.5386 0.6201 0.6337 0.6873 0.7315 0.7443 0.8931 Figure 29: First Modal Shape, a) refers to the classic algorithm: CA, b) refers to the substructured algorithm: CB. ETSEIAT. Fabi´an Lajas Contreras 37 4 APPLICATION FOR A SHELL GEOMETRY Figure 30: Second Modal Shape, a) refers to the classic algorithm: CA, b) refers to the substructured algorithm: CB. Figure 31: Third Modal Shape, a) refers to the classic algorithm: CA, b) refers to the substructured algorithm: CB. ETSEIAT. Fabi´an Lajas Contreras 38 Figure 32: Fourth Modal Shape, a) refers to the classic algorithm: CA, b) refers to the substructured algorithm: CB. 4.5 Discussion of results Inspection of the results shown in the foregoing reveals that, as expected, the natural frequencies computed by the CA and the CB analyses are quite similar (errors below 1 % are obtained for the first six modes). Concerning the deformed shapes associated to each mode, the only significant difference is detected in the third mode, see Figure 13 (wing tail). This discrepancy can be attributed to a simple decompression error, among others. Unfortunately, our analysis fails to reveal the actual reason for such a discrepancy. 5 Application for a Complete Aircraft Configuration 5.1 Chosing the Geometry The aim of this section is to apply the methodology explained in the foregoing section (in which KRATOS and MATLAB are interconnected to automatically generate the required natural frequencies and associated modes) to study the vibration behavior of an entire aircraft structure. To minimize meshing problems, the geometry has been streamlined in several places of the airframe, as explained later in Section 5.3.1 . In the following, we describe the criterion used to guide the choice of the final geometry-. •The geometry should be as close as possible to a conventional commercial airplane, because the specifications are easily accessible —in contrast to military aircrafts . •The geometric model should contain the basic elements present in a conventional airplane. ETSEIAT. Fabi´an Lajas Contreras 39 5 APPLICATION FOR A COMPLETE AIRCRAFT CONFIGURATION Excessively detailed elements such as the landing gear, the engines, lights, control surfaces and pneumatic system will be simply neglected. •The geometry must have two differentiated parts: the inner structure (frames, stringers, spars, ribs ...) and any kind of structural element present in a conventional aircraft. The shell will be formed by the fuselage, the wing tail the stabilizer, and the engine fairings. 5.1.1 Kinds of Models Several alternatives were considered when the project started. More than 10 geometries were selected in order to perform the study. In the following section, the most suitable models will be exposed. The selected geometries can be divided basically into three main kinds of CAD files: •The first kind of geometries consist in classic propeller airplanes. Although there is a huge amount of models available in the Internet, it has to be considered that these geometries tend to contain a huge amount of unnecessary details. Also this kind of airplanes tend to be old military models. Figure 33: First kind of airplanes, classic models A), side view B), profile view C) top view •The second kind of models refer to modern militar jet models. As the previous kind of airplanes, there is also a lot of available material. Nevertheless, this kind of aircrafts have not a trivial structure compared with commercial models. In fact, in this kind of designs, the internal structure is one of the most secret parts. Figure 34: Second kind of airplanes, militar models A), side view B), profile view C) top view •Finally the third kind of models consist in commercial airplanes. These are the preferred models for this project, because as mentioned above, specifications are easier to obtain. Figure 35: Third kind of airplanes, private aviation A), side view B), profile view C) top view ETSEIAT. Fabi´an Lajas Contreras 40 5.1.2 Chosen model Following the guidelines provided in Section 5.1.1 , we choose the third kind of airplane. Specifically, the selected model is a Cessna Citation, a jet commonly used for private aviation. The specifications for this model are shown in Table 7 and correspond to a Cessna Citation XLS model, which is the newest version of this airplane. This updated version contains new details such as different wing configuration, and minor modifications in the rear fuselage. This data have been extracted from [6]. In Figure 36, a representation of the chosen airplane can be seen. Figure 36: Cessna Citation. More information of this product can be found in [5]. Table 7: Cessna Citations XLS Specifications, data courtesy of [6] Exterior Dimensions Wingspan 56 ft 4 in 17.17 m Length 52 ft 6 in 16.00 m Height 17 ft 2 in 5.23 m Interior Dimensions Cabin Height 68 in 1.73 m Cabin Width 66 in 1.68 m Cabin Length 18 ft 6 in 5.64 m Baggage Capacity 800 lb 362.9 kg Weights Maximum Takeoff 20200 lb 9163 kg Basic Operating Weight 12860 lb 5833 kg Useful Load 7540 lb 3420 kg Performance Takeoff Field Length (MTOW) 3560 ft 1085 m Time to Climb FL 450 in 29 min Max Cruise Speed 441 ktas 817 km/h Max Range (Ferry, LRC) 2100 nm 3889 km ETSEIAT. Fabi´an Lajas Contreras 41 5 APPLICATION FOR A COMPLETE AIRCRAFT CONFIGURATION The final design is shown in Figure 47. Material properties are summarized in Table 8. Notice that, given the lack of information concerning the exact geometry of the airplane —and also because of simplicity reasons—, the thickness of all components is considered equal. The distinct substructures employed in modal analysis are sketched in Figure 49. Table 8: Material properties. E[MPa]ν t [mm]ρ[kg/m3] 185 0.29 12 2795.6704 Figure 48: Dirichlet boundary conditions. Figure 49: Structure partition. 5.4 Results As in the previous examples, the number of inner DOFs have been reduced to 6. In table 9, the natural frequencies for each defined substructure are presented. On the other hand, the 3 lowest modal shapes of each substructure are presented in Figure 50, Figure 51 and Figure 52, and the 10 lowest natural frequencies for the total analysis are shown in Table 10. Finally, the 4 lowest modal shapes are displayed in Figure 53, Figure 54, Figure 55 and Figure 56. ETSEIAT. Fabi´an Lajas Contreras 48 Figure 50: Modal Shapes for the first, and the second substructures. ETSEIAT. Fabi´an Lajas Contreras 49 5 APPLICATION FOR A COMPLETE AIRCRAFT CONFIGURATION Figure 51: Modal Shapes for the third and fourth substructures. ETSEIAT. Fabi´an Lajas Contreras 50 Figure 52: Modal Shapes for the fifth and sixth substructures. ETSEIAT. Fabi´an Lajas Contreras 51 5 APPLICATION FOR A COMPLETE AIRCRAFT CONFIGURATION Table 9: Natural frequencies for each substructure Freq [10−2Hz]ω1ω2ω3ω4ω5ω6 sub 1 3.2843 3.2894 3.4046 3.6208 3.8247 4.2056 sub 2 1.7606 4.4526 4.4903 6.4212 6.5100 6.9598 sub 3 1.7506 4.4178 4.4790 6.4196 6.5049 6.9629 sub 4 4.2693 4.5386 4.6705 7.9333 8.9175 11.4050 sub 5 5.9430 6.3922 10.9259 13.2915 14.0920 14.5158 sub 6 5.8864 6.3430 11.1065 13.0727 13.8760 14.6839 Table 10: Natural frequencies. CA refers to Classic Algorithm, whereas CB refers to the CraigBampton algorithm. Freq [10−2Hz]ω1ω2ω3ω4ω5ω6ω7ω8ω9ω10 CA 0.2212 1.1885 1.4333 1.6685 1.7732 1.8537 2.7734 2.8878 3.3185 3.3326 CB 0.2212 1.1870 1.4284 1.6682 1.7709 1.8507 2.7607 2.8739 3.3122 3.3258 Figure 53: First Modal Shape ETSEIAT. Fabi´an Lajas Contreras 52 Figure 54: Second Modal Shape Figure 55: Third Modal Shape ETSEIAT. Fabi´an Lajas Contreras 53 6 CONCLUDING REMARKS Figure 56: Fourth Modal Shape 5.5 Discussion Similarly to the cases studied in previous sections, the agreement between the natural frequencies predicted by the complete and partitioned model is quite acceptable (errors are below 1 %, see table 10), a fact that provides definite evidence that the implementation of the Craig-Bampton algorithm is correct. Concerning the deformed shapes associated to each vibration mode, the resemblance is also rather reasonable —notice that sometimes the vibration modes exhibit the same “shape” yet opposite sign (Figure 56). 6 Concluding remarks •From the values of natural frequencies obtained in three sets of simulations, it may be concluded that the implementation of the substructuring modal analysis methodology —the primary goal of this project— is correct, for the errors in predicting the natural frequencies are relatively low (less than 1 %) and in accordance to values reported in the related literature. Admittedly, discrepancies have been detected in some cases when visually comparing the deformed shapes computed with the classical modal analysis algorithm and the Craig-Bampton method. The exact reasons behind these discrepancies have not been, unfortunately, uncovered by our analyses. Future research should focus on trying to elucidate possible reasons behind these differences. As suggested in the first example, perhaps the causes may lie in the reconstruction process that allows one to plot the deformed shapes. •All simulations carried out have been performed automatically by launching a single script in Matlab. For the shell problems (examples 2 and 3), this script calls a Python program, that in turns, calls the KRATOS software—which is programmed in C++. This software writes the finite element information into binary files that are then read by Matlab again to perform the required modal analysis. Finally, the deformed shapes are plotted by using ETSEIAT. Fabi´an Lajas Contreras 54 GID’s postprocess facilities. In the substructuring version, this task is repeated for each partition. •It should be noted that the simulation of the complete airplane model (design of the geometry, the partitioning, and the final meshing) has required an amount of work and time way higher than the devoted into coding and testing the substructuring algorithm. The geometry developed in this project has been extracted from a private airliner geometry, in which basic aeronautical structures have been added, such as spars, wings and even the torsion box. Even so, the quality of the geometry and meshing cannot be regarded as optimal from a engineering point of view. However, considering the amount of time assigned for this project, it can be considered an acceptable approximation in order to compare the results between the two algorithms studied in this project. 7 Future Lines of Research This project can be the starting point for future research lines, namely •As a first attempt to improve the performance of the program, it could be worthy to check the totality of the code and see if there are any bugs —mainly in the reconstruction process. •Another possible research line would be to optimize the developed code. First of all, some Matlab TM operations are still programmed without vectorization —which is the recommended programming style in matlab. Further improvement can be achieved by parallelizing the modal vibration of all substructures. Another alternative could be simply to translate the entire problem into a general purpose language such as C++ or FORTRAN . This final solution would be, of course, a drastic one. It has to be also considered that, although a general purpose language tend to be faster than Matlab TM, these kinds of languages lack the preimplemented functions that are commonly used in that platform, and this will suppose an extra amount of time in order to implement these auxiliar functionalities. •Of course, it is also possible to simulate a more complex geometry than the ones treated in this study. This kind of models could be obtained by simply creating a finer mesh, as it can be seen in Figure 57, which would require more elements and so more computational effort. The other option is to develop, using CAD technologies a more complex geometry which would contain even more detailed structural components. If a better geometry is desired, it is possible to follow one of the suggested paths or even both, depending on the available resources for the researcher. Figure 57: Mesh refinement: a) 17005 nodes with 36651 elements, b) 35599 nodes with 75739 elements, c) 128794 nodes with 268336 elements. •Trying to increase the flexibility of the code is also an interesting research line. It would consist in changing the code structure at such a way that the geometry of each substructure ETSEIAT. Fabi´an Lajas Contreras 55 REFERENCES may be exchangeable, allowing the simulation of models with different iterations for each substructure, as it is common to see in professional simulations, such the ones performed in the aeronautical industry. References [1] C Felippa. Introduction to finite element methods (asen 5007) course material: Lecture 10, 2013. [2] Eugenio O˜nate. Structural analysis with the finite element method. Linear statics: volume 2: beams, plates and shells. Springer Science & Business Media, 2013. [3] C Felippa. Introduction to finite element methods (asen 5007) course material: Lecture 21, 2013. [4] What is kratos? http://kratos-wiki.cimne.upc.edu/index.php/What_is_Kratos, note = Accessed: 2016-12-25. [5] Cessna citation 2. https://grabcad.com/library/cessna-citation-2-1. Accessed: 20171-16. [6] Cessna citation xls, brochure. http://cessna.txtav.com/~/media/cessna/files/ citation/xlsplus/xlsplus_brochure.ashx, 2015. Accessed: 2016-9-10. [7] Michel G´eradin and Daniel J Rixen. Mechanical vibrations: theory and application to structural dynamics. John Wiley & Sons, 2014. [8] Maurice Petyt. Introduction to finite element vibration analysis. Cambridge university press, 2010. [9] Olgierd Cecil Zienkiewicz and Robert Leroy Taylor. The finite element method: the basis, volume 1. Butterworth-heinemann, 2000. [10] Olgierd Cecil Zienkiewicz and Robert Leroy Taylor. The finite element method: solid mechanics, volume 2. Butterworth-heinemann, 2000. [11] C Felippa. Introduction to finite element methods (asen 5007) course material: Lecture 16, 2013. [12] C Felippa. Introduction to finite element methods (asen 5007) course material: Lecture 17, 2013. [13] C Felippa. Introduction to finite element methods (asen 5007) course material: Lecture 18, 2013. [14] Yair M Altman. Accelerating MATLAB Performance: 1001 tips to speed up MATLAB programs. CRC Press, 2014. [15] Design breakdown of the p-51 mustang. http://migrate.legendsintheirowntime.com/ LiTOT/P51/P51_IA_4407_DA.html. Accessed: 2016-12-29. [16] The engineering toolbox. http://www.engineeringtoolbox.com/poissons-ratio-d_1224. html. Accessed: 2016-12-29. ETSEIAT. Fabi´an Lajas Contreras 56 [17] Michael Bauccio et al. ASM metals reference book. ASM international, 1993. [18] C Esbr´ı. A-320: Structures, 2015. [19] Eduardo Pinheiro, Ricardo Bianchini, Enrique V Carrera, and Taliver Heath. Load balancing and unbalancing for power and performance in cluster-based systems. In Workshop on compilers and operating systems for low power, volume 180, pages 182–195. Barcelona, Spain, 2001. A Costs A.1 Workings Hours To calculate the cost of time, an hypothetical salary of 36000 AC has been considered. Using a working week of 40 hours, total hour-rate would be 18.75 AC/h. Starting by this assumption and the hours listed below, the total amount corresponding to working hours can be seen in Table 11. Table 11: Working hours costs Concept Time Cost Information research 90 1687.5 Development of the Project Charter 10 187.5 CAD software development 333 6243.75 GID file reverse engineering 15 281.25 Kratos files reverse engineering 20 375 Generation of the first problem geometry 10 187.5 Generation of the second problem geometry 24 450 Generation of the third problem geometry 264 4950 Data preparation and postprocess 40 750 Generate the mesh 5 93.75 Prepare external geometry to import 5 93.75 Modify Kratos source to extract matrices 30 562.5 Program implementation 214 4012.5 Implement vectorized modal analysis 20 375 Implement Craig-Bampton code 100 1875 Create Kratos interface 16 300 Implement beam element 20 375 Implement GiD geometry creator 48 900 Implement Kratos reader 10 187.5 Final data analysis and postprocess 15 281.25 Extract the results for the first problem 10 187.5 Compare the results of both problems 5 93.75 Documents and others 50 937.5 TOTAL 752 hours 14100 AC ETSEIAT. Fabi´an Lajas Contreras 57