- Email: [email protected]
Optimization approach for identification of dynamic parameters of localized joints of aircraft assembled structures
JID:AESCTE AID:4121 /FLA
[m5G; v1.221; Prn:27/07/2017; 8:46] P.1 (1-12)
Aerospace Science and Technology ••• (••••) •••–•••
1
67
Contents lists available at ScienceDirect
2
68
3
69
Aerospace Science and Technology
4
70
5
71
6
72
www.elsevier.com/locate/aescte
7
73
8
74
9
75
10
76
11 12 13
Optimization approach for identification of dynamic parameters of localized joints of aircraft assembled structures
14 15
77 78 79 80
a,∗
b
a
a
a
a
81
17
S. Hernández , E. Menga , S. Moledo , L.E. Romera , A. Baldomir , C. López , M. Cid Montoya a
18
a
84
16
19
Structural Mechanics Group, School of Civil Engineering, Universidade da Coruña, A Coruña, Spain b Component Loads and Aeroelastic Department – EGLRG, Airbus Operations Getafe, Getafe, Spain
20
25 26 27
i n f o
a b s t r a c t
88 89
Article history: Received 12 September 2016 Received in revised form 25 May 2017 Accepted 17 July 2017 Available online xxxx
28 29 30 31 32 33
85 87
a r t i c l e
23 24
83
86
21 22
82
Keywords: Aircraft design Assembled structures Ground vibration test Dynamic analysis Optimization methods
34 35 36 37 38 39 40 41 42
Ground Vibration Test (GVT) is a crucial step in the design process of a new aircraft. The test involves the analysis of different load configurations that helps to understand the real structural behaviour, to certificate some extreme load cases and also to validate and improve the Finite Element (FE) model used in the design. The validation of the FE model is done by making a comparison between Frequency Response Functions (FRF) from the GVT and the FE model at some response points. Ideally, numerical and experimental results should be in agreement. This is certainly true in the case of single components but in the case of assemblies it is common that numerical and experimental results show discrepancies, especially in the case of dynamic analysis. These differences increase as the number of components in the assembly grows, because connections are usually not represented sufficiently well in the FE model. The components of these assembled structures are usually connected at a small number of points, and they can be seen as local sources of stiffness and damping. In the FE model these joints can be represented as lumped spring/damper elements in a simplified way. In this paper a formulation of the dynamic equilibrium equation employing a decoupled approach for the joint mechanical properties of an assembly is presented. Later, an optimization process based on this formulation that involves the joint properties as design variables has been developed in order to identify the joint parameters that improve the adjustment between numerical and experimental FRF in the case of harmonic excitations. The approach has been applied to a FE model that represents the connection between a segment of fuselage and the tail cone and to a real aircraft system where experimental data is available. The results achieved with this procedure are very efficient and promising. © 2017 Elsevier Masson SAS. All rights reserved.
90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108
43
109
44
110
45 46
111
1. Introduction
47 48 49 50 51 52 53 54 55 56 57 58 59
In the manufacturing process of an aircraft its components are built individually, even at different locations. Wings, fuselage sections, vertical or horizontal stabilizers and other components are assembled after their production in a final assembly line. These components are usually connected at few locations to create an assembly. Fig. 1a shows an example of a single component of an aircraft fuselage while Fig. 1b presents an assembly of two components. In the case of single components, modelling and dynamic prediction has been well developed, and accurate predictions may be achieved. However, in the case of assemblies the quality of the results decreases with the complexity of the connection. This arises
60 61 62 63 64 65 66
*
Corresponding author. E-mail address: [email protected] (S. Hernández).
http://dx.doi.org/10.1016/j.ast.2017.07.026 1270-9638/© 2017 Elsevier Masson SAS. All rights reserved.
due to two possibilities: the component model is not sufficiently accurate or the connection has not been represented satisfactorily. The second reason is frequently the most probable. The appropriate characterization of mechanical properties in assemblies is a recurrent problem for both static and dynamic analyses. Fitting the correct value of stiffness and local damping involved in this type of connections is crucial to obtain accurate models able to predict reliable responses to any kind of loads. To efficiently analyse the dynamic behaviour of assembled structures several methods using substructure coupling have been developed. The Component Mode System (CMS) method was proposed by Hurty [1], and has been widely used and developed [2–4] for the last forty years to predict the frequencies and modal behaviour of complex FE models on the basis of the subsystems modal properties that take part of them. In the case of using experimental Frequency Response Functions (FRF) as response data, the FRF based coupling method proposed by Jetmundsen [5] en-
112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE
2
AID:4121 /FLA
[m5G; v1.221; Prn:27/07/2017; 8:46] P.2 (1-12)
S. Hernández et al. / Aerospace Science and Technology ••• (••••) •••–•••
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39
Fig. 1. Example of components and assembled structures.
40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
ables the structural synthesis in the frequency domain. Both methods allow combining experimental and numerical results and were initially developed considering rigid joints connections between assembled parts. In order to improve the dynamic analysis of assembled structures considering more realistic joint connections, several modifications of these two methods have been proposed, like the General Joint Description Method (GJDM), to take into account the effect of elastic media between subsystems [6,7]. Other approach to study the dynamic response of assembled structures is the Structural Dynamic Modification Method (SDMM), which allows evaluating the incidence of a set of modifications in the dynamic response of the structure without the need to perform FE simulations repeatedly. The SDMM is particularly efficient when the modifications are concentrated in lumped elements. The modified structural response can be expressed as a function of the baseline FE model dynamic database and the set of modifications. A good theoretical approach of SDMM and some work based on it is presented in [8–10]. From a modelling point of view connections can be seen as local sources of stiffness and damping, in most cases with a nonlinear response. Connections are usually performed with prestressed bolts, rivets, pins, adhesives or elastomeric parts, and the complex behaviour of these elements with contact forces, friction, precompression loads and nonlinear material laws substantially modifies the dynamic response of the global system due to changes in stiffness and damping [11–13]. The accurate characterization and mod-
elling of the dynamic behaviour of these connections is a prolific field of research and important efforts have been made in both the numerical modelling and parameter identification by comparing with experimental results [14–17]. They are usually defined by nonlinear constitutive laws, which parameters might depend from a great variety of aspects, such as load application points, load frequency or magnitude, relative displacement value or load duration. Detailed three-dimensional joint models are accurate but not applicable in the analysis of large structures with a high number of connections [18]. The most commonly used models are lumped parameters with the joint modelled as a single point [19–21] and to a lesser extent thin layer elements [22]. Obtaining nonlinear dynamic responses is a difficult task, especially for large structures, due to the necessity of employing iterative methods in time or frequency domain [23] and the complexity of nonlinear vibrations. Among the solutions currently employed to solve this problem, some can be highlighted such as the Harmonic Balance Method (HBM) in the frequency domain [22,24]. HBM is a general method based on identifying periodic solutions as truncated Fourier series with a variable number of harmonics; simplified models for nonlinear time history responses [25], or nonlinear internal force vectors expressed as an equivalent stiffness matrix with dependency of vibration amplitude [26]. During the project of a new aircraft a Finite Element (FE) model is created based on expert knowledge, FE modelling state-of-theart and information learnt from previous successful designs. Once all parts of the new aircraft are designed and manufactured a Ground Vibration Test (GVT) is performed. The GVT is done for validating and improving the structural model. During the GVT several shakers or Driving Points (DP) are located at specific positions of the structure in order to excite the full aircraft, and a great amount of measurement devices, like accelerometers or Response Points (RP) are located throughout the whole aircraft to register as much information as possible. Fig. 2 shows a GVT carried out in a civil aircraft and in a military aircraft. The GVT is carried out using different load configurations. Meanwhile the GVT is performed the FE model is updated. What is learnt from the experiment is introduced in the model, and predictions from the model are used to obtain better experimental results. Since the test is performed in the final stages after the assembly of the prototype, one of the priorities is to perform testing and model updating in the shortest period of time. Updating the mathematical model is crucial to make reliable flutter predictions for flight tests [27–29], and this is why recognizing experimentally the low-frequency modes of the whole aircraft structure is the first scope [30]. It is also an objective of GVT the certification of Fan-Blade-Off and Windmilling analysis [31], because these extreme cases can never be tested under real conditions. Additional testing is required in the case of complex aircraft designs trying to quantify the non-linear behaviour of the structure comparing the obtained responses with different load levels. In the last years new issues have arisen like passenger’s comfort. In the past the use of the Phase Resonance Method or socalled Normal Mode Testing has been used for GVT of large aircrafts, which essentially consist of applying single sine excitations at the structural natural frequencies. This method is quite good for separating closely spaced modes, but has a very time-consuming testing procedure. Hence, it has been substituted by the so-called phase separation techniques that find the aircraft modes by evaluating the FRFs. The signals used in this method are harmonic, periodic signals like multi-sines, swept-sines [32] or periodic chirp, transient signals like impulses or random. The main target of matching FE models and GVT is to validate the FE model in order to predict the behaviour of the structure. However it is not uncommon to find many discrepancies between both results. For many years the lack of reliability when modelling
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE AID:4121 /FLA
[m5G; v1.221; Prn:27/07/2017; 8:46] P.3 (1-12)
S. Hernández et al. / Aerospace Science and Technology ••• (••••) •••–•••
3
where m, c, and k are the mass, modal damping and stiffness matrix, while u and f represent the displacement and the external forces vectors. In the case of assembled structures, the stiffness matrix represents the stiffness of the whole assembly. Since this work pursues to improve the fitting of the FE results against experimental testing by obtaining the adequate values of the mechanical properties of the joints that connect the components of the assembly, it would be necessary to isolate such properties from the dynamic equilibrium equation. In this sense, a novel formulation with the aim of identifying the dynamic parameters of the joints in assembled structures in order to fit numerical FRF with experimental FRF is developed in this section for being applied later on. As seen in Hernández [33], assuming linear spring joints, the stiffness matrix k can be divided in two parts (ks + kns ) based on the SDMM:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 18 19 20 21 22 24 25 26 27 28 29 31
(2)
where the matrix k g has the products of the structural damping by the variable spring stiffness:
32 33
kg =
35
⎛
Fig. 2. Ground Vibration Test of a civil and a military aircraft.
37
41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59
the junctions in complex assemblies has been underestimated. This research presents a decoupled formulation of the dynamic equilibrium equation that isolates the local stiffness and damping of the FE model connections from the rest of the structure, aiming to identify through an optimization process their optimum values in order to improve the performance of FE models of assembled structures for reaching a good agreement between experimental and numerical results. The novelty of this research lies in the formulation exposed in Section 2, which was developed in the framework of an industrial collaboration with Airbus Operations for improving the correlation between the experimental and numerical results of assembled structures connected at a small number of joints. In this sense, the dynamic behaviour of the parts of the assembly is assumed to be well represented in the FE analyses, being the focus of the discrepancies in the numerical behaviour of the joints, which is certainly the most usual case in real situations. This methodology is successfully applied to a simplified FE model that represents the assembly between a rear fuselage section and the tail cone of an aircraft and to an aircraft system that consist of a Vapour Cycle Refrigeration Unit (VCRU) connected to the fuselage of an aircraft. The latter case is studied considering real experimental data.
60 61 62 63 64
2. Formulation of the dynamic analysis of assembled structures through a decoupled approach The traditional formulation of the dynamic equilibrium of a structure is:
65 66
71 72 73 74 75 76 77 78 79 80 81 83 84 85 86 87 89 90 91 92 93 94 96 97 98 99
34
40
70
95
k = (kns + ks )(1 + ig ) = kns + ks + ik g
30
39
69
88
In a limit case, the entire spring stiffness might be assigned to ks , while kns would contain just the stiffness of the separate components. However, in practical applications this would worsen the behaviour of the assembly since the components will not even be partially connected. Considering structural damping g, the total structural stiffness can be split into the sum of three matrices.
23
38
68
82
– A matrix kns that contains the stiffness of all the separate parts that builds the assembly and a partial and constant part of the spring stiffness. – A matrix ks than contains the rest of the uncoupled spring stiffness, considering it as variable.
17
36
67
¨ + cu˙ + ku = f mu
(1)
100
.. ⎜ . ⎜. . . (k + k ) g nsi si i ⎜ ⎜ −(k + k ) g nsi si i ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ . . .
101
⎞
. . .
⎟ ⎟ ... p −(knsi + ksi ) g i . . . ⎟ ⎟ ...q (knsi + ksi ) g i ⎟ ⎟ . . ⎟ . . .. . ⎟ . . ⎟ . . . (knsj + ksj ) g j −(knsj + ksj ) g j . . . ⎟ ⎟ ...r ⎟ ...s −(knsj + ksj ) g j (knsj + ksj ) g j ⎟ ⎠ . .. . . . . . . . . . . . .
p
q
r
s
(3) In which p and q, r and s, are the degrees of freedom that connect the springs i and j respectively, as seen in Equation (3). The total stiffness of the joint, represented as ki is the sum of k si and knsi , and g i is the structural damping value of the joint, which is the local damping parameter of the joint. The matrix ks has the variable stiffness information of the joints:
⎛.
⎜. . . k −ksi ⎜ si ⎜ −ksi k si ⎜ ⎜ .. ⎜ ks = ⎜ . ⎜ ⎜ ⎜ ⎝ .. .
.. .
p
q
... .. . . −ksj . . . ksj −ksj k sj .. . .. .. . . ..
r
s
⎟ ... p ⎟ ⎟ ...q ⎟ ⎟ ⎟ ⎟ ⎟ . . .⎟ . . .r ⎟ ⎠ ...s .. .
103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
⎞
.. .
..
102
121 122 123 124
(4)
125 126 127 128 129 130 131 132
JID:AESCTE
AID:4121 /FLA
[m5G; v1.221; Prn:27/07/2017; 8:46] P.4 (1-12)
S. Hernández et al. / Aerospace Science and Technology ••• (••••) •••–•••
4
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
where p and q are the joints that define the i spring. Thus, the dynamic equilibrium under a vector of forces f can be rewritten as:
20 21 22
Defining the displacement vector u as a function of the spectral matrix which contains the natural vibration modes, u becomes:
u = q,
(6)
where is the spectral matrix that contains the eigenvectors associated with the matrix kns normalized to the total mass. The vectors f and q can be written as
q = ve i ωt .
(7)
Introducing these expressions in (5) the dynamic equilibrium becomes
Iq¨ + Cq˙ + Kns q + Ks q + iK g q = F,
23
C = Φ cΦ
24
and
25 26 27 28 29
32
T
F = Φ p,
Kns = Φ T kns Φ
⎛
35 36 37 38 39 40 41 42 43 44 45 46 47 48
Ks = Φ T ks Φ
⎜
ω12 . . .
. Kns = ⎝ ..
..
K g = Φ T k g Φ,
⎞
⎛
.. ⎟ . ⎠
C=⎝
⎜
. 2 . . . ωm
33 34
(9)
(10)
where the matrices C and Kns are easily built from the modal information as follows,
30 31
(8)
where T
2ζ1 ω1
.. .
⎞
... ⎟ .. .. ⎠ , (11) . . . . . 2ζm ωm
being m the number of natural frequencies. The vectors p, F and v in general are complex, so
p = pr + ipi
F = Fr + iFi
v = vr + ivi .
Therefore, it can be written that
Rns + Ks ωC + K g
−(ωC + K g ) Rns + Ks
vr vi
=
Fr Fi
(12)
(13)
,
where
Rns = Kns − ω2 I.
(14)
Solving the system of linear equations (13) the vector of dis¨ can be obtained placements u and the vector of accelerations u for each type of dynamic load as
49 50 51 52 53 54 55 56 57 58 59 60 61
u = q = |vr + ivi |e
64 65 66
¨ = −ω |vr + ivi |e u 2
i ωt
.
(15)
A specific displacement uk or acceleration u¨ k will be written as
uk = ψ k vr + i (ψ vi ) e i ωt
u¨ k = −ω ψ k vr + i (ψ vi ) e i ωt , (16) 2
where ψ k is the k-th row of matrix . This can be expressed as,
uk = (ψ k vr )2 + (ψ k vi )2 e i (ω+θk )t = Ae i (ω+θk )t u¨ k = −ω2 (ψ k vr )2 + (ψ k vi )2 e i (ω+θk )t = −ω2 Ae i (ω+θk )t
(17) (18)
where θk is the phase angle with a value of
62 63
i ωt
θk = arctg
ψ k vi ψ k vr
69
(5)
18 19
68 70
¨ + cu˙ + ku = mu¨ + cu˙ + (kns + ks + iks )u = f. mu
f = pe i ωt
67
(19)
and therefore A is the magnitude of the FRF of displacement uk for the load with frequency ω .
71
Fig. 3. Springs definition.
This formulation allows to obtain the displacements, velocities and accelerations of a structure in a novel way. It is only required the spectral matrix ns , without the variable springs stiffness, and the natural vibration frequencies. The main advantage is its flexibility. From the global modal information, which only needs to be obtained once, several fast trials with different springs definition can be made by just changing the matrices k g and ks . The analysis of one FRF can be done following these steps: – Analyse the eigenvalue problem and build the spectral matrix ns , after removing the variable spring stiffness. – Assemble of stiffness matrices k g and ks . – Compute K g and Ks as seen in expression (10). – Create the matrix Rns as seen in expression (14). – Solve the system of equations in (13), and compute expression (16) to obtain the displacements and accelerations. The final purpose of this formulation is to produce a fitting between numerical and experimental results. The main advantage is the fact that the truncated spectral matrix ns remains constant, which simplifies the process, removing the need of solving the dynamic modal analysis of the full structure every time a modification in the joints is made. These modifications affect only the matrices k g and ks .
72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98
3. Validation of the uncoupled formulation
99
The formulation proposed to solve this problem divides the stiffness matrix into the sum of three as seen in expression (2). In order to validate this formulation several numerical tests with different ratio of uncoupled springs stiffness have been made. The analysis has been done in the structural laminar model defined in Fig. 3, which represents the connection between two sections of an aircraft. The model consists of an 8 meters long cylinder of 4 meters of radius and a truncated cone of 4 meters long and varying radius from 4 to 2 meters. The cylinder that represents the fuselage is considered built-in and a sinusoidal unitary load is applied at the tip of the cone that represents the tail of the fuselage, as shown in Fig. 4. A constant value of 2 mm has been applied to the shell elements of the model and aluminium, with Young’s modulus value of 70000 MPa and Poisson’s ratio of 0.2, has been selected. The connection between the cylinder and the cone is composed by 5 joints, with stiffness in 9 degree of freedom in total. The position of the joints can be identified by the angle α as presented in Fig. 5, which takes the following values: 45◦ , 135◦ , 225◦ , 270◦ and 315◦ . The location and connectivity of the joints can be seen as well in Fig. 5. These joints represent the assembly between a segment of fuselage and the tail cone of a real aircraft. In order to validate the formulation several analyses with different coupled stiffness ratio have been carried out for coupled stiffness in range 40% to 100%. For instance, if a ratio of 60% of coupled stiffness is defined it means that 40% of the stiffness of the lumped elements is in the matrix ks (matrix of variable springs stiffness). The remaining stiffness is obtained from the modal analysis of the finite element model changing the stiffness of the joint spring elements to the 60% of their initial value. This stiffness is computed directly using the matrix Kns with no more information than the eigenvalues as seen in (11).
101
100 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE AID:4121 /FLA
[m5G; v1.221; Prn:27/07/2017; 8:46] P.5 (1-12)
S. Hernández et al. / Aerospace Science and Technology ••• (••••) •••–•••
5
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11 12
77
Fig. 4. Model geometry and load location.
78
13
79
14
80
15
81
16
82
17
83
18
Fig. 7. Detail of the comparison of the uncoupled formulation at node 5432 DOF 2, in the frequency range from 15 to 19 Hz, without considering structural damping.
19 20
The formulation shows accurate results at different ratios of coupled stiffness when compared with the 100% coupled solution obtained in MATLAB. Moreover, a validation of the formulation was accomplished after comparing the results of the 100% coupled solution with the numerical one obtained with MSC NASTRAN through a SOL111, a modal frequency response analysis. In Fig. 7, which is a detail of Fig. 6 in the frequency range from 15 to 19 Hz, the effects of modifying the ratio of coupled stiffness can be recognized. It is clear that with higher values of coupled stiffness better results will be reached, being the 100% case of coupled stiffness the exact solution. Fig. 8 and Fig. 9 are the same as Fig. 6 and Fig. 7 but considering added structural damping with a value of 0.1. Once again the comparison with MSC NASTRAN has been done, and a total superposition was obtained. Even at low levels of coupled stiffness the formulation shows to be quite precise.
22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
Fig. 5. Location and connectivity of joints.
38 39 40 41 42 43 44 45 46 47 48 49 50
85 86
21
37
84
The structure has a constant modal damping of 2%, whereas each spring has a stiffness value of 5000 N/m. The damping ratio of 2% is based on experimental vibration studies found in literature [34,35]. Modal analysis has been performed with MSC NASTRAN [36] to obtain the spectral matrix and the natural frequencies on a range from 0 to 100 Hz, which included 229 modes. The formulation developed in the previous section has been implemented in MATLAB [37]. The results shown in Figs. 6 to 9 were obtained with a vertical and unitary harmonic force located at node 5432 (Fig. 4) and the response studied was the vertical acceleration at the same node.
87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103
4. Identification of mechanical parameters by an optimization process
104
In assembled structures certain parameters such as the stiffness or the local damping related to joints may be unknown. In this research the fitting between experimental results obtained from the GVT and numerical results from the FE model has been done by identifying these parameters. Thanks to the uncoupled formulation any change can be easily made in the stiffness matrix without any further information from the FE model. Just varying the uncoupled matrix ks a better fitting between experimental information and numerical results can be achieved.
107
105 106 108 109 110 111 112 113 114 115 116
51
117
52
118
53
119
54
120
55
121
56
122
57
123
58
124
59
125
60
126
61
127
62
128
63
129
64
130
65 66
131
Fig. 6. Comparison of the uncoupled formulation at node 5432 DOF 2 without considering structural damping.
132
JID:AESCTE
AID:4121 /FLA
[m5G; v1.221; Prn:27/07/2017; 8:46] P.6 (1-12)
S. Hernández et al. / Aerospace Science and Technology ••• (••••) •••–•••
6
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
Fig. 8. Comparison of the uncoupled formulation at node 5432 DOF 2 considering a structural damping of 0.1.
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27 29 30 31 32 33 34 35
Fig. 9. Detail of the uncoupled formulation at node 5432 DOF 2, in the frequency range from 15 to 19 Hz, considering a structural damping of 0.1.
36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
93
Fig. 10. Strategy of the optimization process.
28
Furthermore, an optimization process can be defined to find the best value of the objective function. The design variables of the problem are the parameters involved in the connection (spring stiffness and structural damping) and also a set of modal damping values that cover the ranges of the frequencies involved. The objective function for this optimization problem was defined as the sum of two terms, one related to the natural frequencies and another related to the FRF value at those frequencies. Minimization of differences between the numerical and experimental natural frequencies is done in the first term of the objective function, approaching horizontally the numerical and experimental FRF. The second term of the objective function aims to fit the FRF vertically by minimizing the difference between the peaks values, which are the FRFs evaluated at their natural frequencies. The objective function (20) is formulated as the sum of squares of the differences between numerical i-th natural frequency ωiN U and experimental i-th natural frequency ωi E X plus the sum of squares of the differences between numerical values of the FRF evaluated at the i-th natural frequency, FRF N U (ωiN U ), and the experimental value of the FRF evaluated at the i-th natural frequency, FRF E X (ωi E X ).
is the i-th experimental natural frequency, and ωiN U is the i-th numerical natural frequency. In expression (20) the experimental values extracted from the GVT do not vary since they are the target values of the optimization process. The remaining terms are evaluated during the optimization process, because any change in a parameter would change the numerical FRF, not only the position of the natural frequencies but also the peak values at those frequencies. Fig. 10 shows how the fitting between numerical and experimental is reached. The numerical natural frequencies are defined as the points where the maximum values of the FRFs are located. To start the optimization process these natural frequencies are obtained from the solution of the eigenvalue problem. Once all numerical natural frequencies are known a range of frequencies ω is defined in order to reduce the analysis time and evaluate the FRF only around the desired frequencies. Initially, ranges are defined considering the numerical and experimental natural frequencies. The optimization process assures that any changes in the design variables will produce numerical results that are included within the interval [min( w iN U − w , w i E X − w ), max( w iN U + w , w i E X + w )]. Expression (20) is formulated for a single FRF. It might be more suitable to fit simultaneously as many FRFs of degrees of freedoms as possible in order to enhance the accuracy between numerical and experimental results. A more general objective function that includes more response points, and consequently more FRFs, can be defined as
F=
n
F= (ωi E X − ωiN U )2
N FRF
+
2 FRF E X (ωi E X ) − FRF N U (ωiN U )
N w
λj (ωi j E X − ωi jN U )2
j =1
i =1 n
+ (20)
i =1
where n is the number of peaks considered in the optimization, FRF E X is the experimental FRF, FRF N U is the numerical FRF, ωi E X
94
i =1
Nw
95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125
2 FRF j E X (ωi j E X ) − FRF jN U (ωi jN U ) ,
126
(21)
i =1
where N FRF is the number of FRF taken into account, λ j is a weight factor that has been set to 1 in this study, N w is the number of frequencies considered, FRF j E X is the j-th experimental FRF, FRF jN U
127 128 129 130 131 132
JID:AESCTE AID:4121 /FLA
[m5G; v1.221; Prn:27/07/2017; 8:46] P.7 (1-12)
S. Hernández et al. / Aerospace Science and Technology ••• (••••) •••–•••
7
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14 15
80
Fig. 11. Fitting of FRF of vertical acceleration at node 5432.
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31 32
97
Fig. 12. Location of the parameters taken into account in the optimization case with 4 FRF response points simultaneously.
33 34 35 36
99
is the j-th numerical FRF, ωi j E X is the i-th experimental natural frequency of the j-th FRF, and ωi jN U is the i-th numerical natural frequency of the j-th FRF.
37 38 39
5. Application example modelling an assembled aircraft structure
40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The formulation exposed previously has been applied to the structural model defined in Fig. 4, which represents the connection of two fuselage segments of an aircraft. Fig. 5 shows the connectivity of this section, composed by 5 springs working in 9 DOF in total. The experimental data has been simulated by changing some parameters of the model, such as the springs stiffness or the modal damping. Then a frequency domain analysis has been carried out using MSC NASTRAN. The results of this analysis are used as the objective ones in the optimization process, simulating the experimental results obtained from a GVT. Two different configurations have been taken into account in this model: one considering as design variables the spring’s stiffness and the modal damping of several natural frequencies and another considering in addition as design variables the local damping modelled as structural damping. The optimization problem was solved using MATLAB with a Sequential Quadratic Programming (SQP) algorithm. In both cases the dynamic excitation was a harmonic unitary force applied in vertical direction at node 5432, on the top tip of the cone (Fig. 4).
61 62 63
5.1. Optimization without including structural damping as design variable
64 65 66
98
Firstly, an optimization problem with only one response point was solved. This problem uses the objective function defined in
expression (20). The optimization process was formulated considering as design variables the stiffness of the 5 springs, which are 9 design variables (one for each DOF), plus 9 modal damping coefficients selected at the natural frequencies excited in this response. The response considered was the FRF of the vertical acceleration at node 5432. The spring stiffness values of the initial design are a 90% of the spring stiffness values employed in the objective design. Moreover, the values of the damping were modified from a 2% in the objective design to a 3% in the initial design. The objective design in this academic example is obtained numerically through NASTRAN, however in practical applications it comes from experimental tests. These values of initial and objective parameters were used in all cases studied in this research. The number of eigenvectors used is the same as in the validation example. The proposed optimization process is able to reduce the objective function from an initial value of 51.8271 to a final value of 0.8668, which means an improvement of 98.32%. The initial, objective and optimized FRF of the response point selected can be seen in Fig. 11. The tolerance established for the objective function was a relative error of 1 × 10−6 in all cases and the algorithm performed 128 iterations before convergence in this example. After seeing the good results achieved by the optimization process a new problem considering 4 different response points was formulated. The final purpose is to use as many response points as possible, which would mean a better accuracy between FE model and experimental results. In this case the selected objective function was the one defined in expression (21). The following RPs were used:
100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128
– Node 5432 DOF 2, which is the vertical direction of the node where the load is applied. – Node 3405 DOF 3, which corresponds to the longitudinal direction of the upper left spring.
129 130 131 132
JID:AESCTE
8
AID:4121 /FLA
[m5G; v1.221; Prn:27/07/2017; 8:46] P.8 (1-12)
S. Hernández et al. / Aerospace Science and Technology ••• (••••) •••–•••
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
Fig. 13. Fitting of vertical acceleration at node 5432 with 4 response points in parallel.
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29 30
95
Fig. 14. Fitting of longitudinal acceleration at node 3405 with 4 response points in parallel.
31 32 33 34 35 36
97
– Node 3443 DOF 1, which corresponds to the horizontal direction of the bottom spring. – Node 3455 DOF 2, which corresponds to the vertical direction of the bottom right spring. The position of the parameters used in this optimization case can be seen in Fig. 12.
37 38 39 40 41 42 43 44 45 46 47 48 49
In the example, whose results are provided from Figs. 13 to 16, 25 design variables have been used; 9 matching the springs stiffness as seen previously in the one parameter optimization case and 16 corresponding to the modal damping coefficients. In this case there are 5 more than in the previous example because some natural frequency might be excited in one parameter and might not be in another one. To achieve good results, all the modal damping coefficients of the modes excited at any response point must be considered. In this problem the results were also very accurate. The initial value of the objective function was 77.1348 and the final value was 5.0214, which means an improvement of 93.49%.
50 51 52
5.2. Optimization including structural damping as design variables
53 54 55 56 57 58 59 60 61 62 63 64 65 66
96
In this example local damping, modelled as structural damping, has been added to the optimization process. Consequently, 9 new design variables were added to the problem, one at each spring component. The initial values of all structural damping are the double of the ones employed in the objective design, which as exposed in section 5.1 represent the experimental data. Following the same idea presented in section 5.1, a simpler optimization considering just one response point was initially formulated. In this case the total number of design variables was 27. The optimization process shows also good results with the structural damping formulation as it can be seen in Fig. 17, in fact the improvement was a 95.05%. On the following pictures g refers to structural damping. The initial point of the optimization process had an ob-
jective function value of 40.2751, and after optimization it was 1.9936. The next step was to carry out the formulation with 4 response points in parallel. This resulted in a bigger optimization problem in which 31 design variables were used. The response points used in this case were the same used in the example of section 5.1. In this example the decrease of the objective function is a bit lower than in previous cases, but as it can be seen in Figs. 18–21 the structural behaviour has drastically been improved. An improvement of 91.97% was achieved, reducing the initial value of the objective function from 59.9372 to 4.8110. It must be pointed out that the response with less magnitude, such as node 3405 and mainly node 3443, have lower improvements that the response points with larger responses (5432 and 3455), maintaining higher discrepancies in comparison with more excited RP. This is due to the use of absolute error in the formulation of the objective function, which was defined intentionally because the main target is to improve the larger responses. However, this can be modified by considering relative instead of absolute errors in the objective function or weighting the term of the objective function associated with the most important responses. In Table 1 there is a summary of the results achieved in this research.
98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121
6. Application to an aircraft system with experimental data available
122 123 124
The methodology presented in Section 4 has been applied to the FE model of an aircraft system for which experimental data is available. The test is carried out to assure an appropriate level of vibrations in the structure. The system is a Vapour Cycle Refrigeration Unit (VCRU), which is usually located in pairs at each side of the fuselage, being connected to it through a suspension system composed by a frame and a set of struts, each one containing a rubber mount that is modelled as a lumped spring
125 126 127 128 129 130 131 132
JID:AESCTE AID:4121 /FLA
[m5G; v1.221; Prn:27/07/2017; 8:46] P.9 (1-12)
S. Hernández et al. / Aerospace Science and Technology ••• (••••) •••–•••
9
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14 15
80
Fig. 15. Fitting of transversal acceleration at node 3443 with 4 response points in parallel.
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31 32
97
Fig. 16. Fitting of vertical acceleration at node 3455 with 4 response points in parallel.
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48 49
114
Fig. 17. Fitting of vertical acceleration at node 5432 with structural damping considering only this response point.
115
50
116
51
117
52
118
53
119
54
120
55
121
56
122
57
123
58
124
59
125
60
126
61
127
62
128
63
129
64
130
65 66
131
Fig. 18. Fitting of vertical acceleration at node 5432 with structural damping and 4 response points in parallel.
132
JID:AESCTE
AID:4121 /FLA
[m5G; v1.221; Prn:27/07/2017; 8:46] P.10 (1-12)
S. Hernández et al. / Aerospace Science and Technology ••• (••••) •••–•••
10
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
Fig. 19. Fitting of longitudinal acceleration at node 3405 with structural damping and 4 response points in parallel.
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
Fig. 20. Fitting of vertical transversal at node 3443 with structural damping and 4 response points in parallel.
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
Fig. 21. Fitting of vertical acceleration at node 3455 with structural damping and 4 response points in parallel.
48 49 50 51
114
Cases
52 53 54
115
Table 1 Summary of all optimizations carried out.
1 RP 4 RP
116
Without structural damping
117
With structural damping
Initial F obj
Final F obj
%
Initial F obj
Final F obj
%
51.8271 77.1348
0.8684 5.0214
98.32 93.49
40.2751 59.9372
1.9936 4.8110
95.05 91.97
55
58 59 60 61 62 63 64 65 66
119 120 121
56 57
118
122
element. The components of the system (compressor, evaporator, fan, . . .) are modelled as concentrated masses that are connected to the frame through rigid elements and springs. Fig. 22 shows the FE model of a VCRU and their location within the fuselage. The response points selected for the experimental validation are the point where the excitation is applied (RP1) and two points of the frame located in the upper side of the frame, near the connection to the fuselage (RP2 and RP3). The locations of the response points are shown in Fig. 23.
The excitation in the experimental test is introduced as a frequency dependent and periodically rotating radial force that is applied at RP1. The accelerations are measured in the three response points at five excitation frequencies. The experimental results provided are a modified version of the results of a real project. In this example, the target is to fit the numerical FRF to the experimental data at the frequency points where information is available. The design variables considered are the modal damping of the structure and the local stiffness and damping of the shock mounts. The modal damping is frequency dependant and changes
123 124 125 126 127 128 129 130 131 132
JID:AESCTE AID:4121 /FLA
[m5G; v1.221; Prn:27/07/2017; 8:46] P.11 (1-12)
S. Hernández et al. / Aerospace Science and Technology ••• (••••) •••–•••
11
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16 17
Fig. 22. FE model of the VCRU and location within the fuselage.
82 83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29 30
95
Fig. 23. Location of the response points selected.
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47 48
Fig. 24. Fitting of the experimental results and the numerical FRF at three RP in parallel. Response in RP1.
Fig. 25. Fitting of the experimental results and the numerical FRF at three RP in parallel. Response in RP2.
49
116
7. Conclusions
51 52
56 57 58 59 60 61 62 63 64 65 66
117 118
53 55
114 115
50
54
113
every 10 Hz in the range of study [0–200] Hz, leading to 21 damping design variables. Figs. 24–26 show the experimental results measured in the three response points at the target excitation frequencies, the initial numerical FRF and the numerical FRF obtained after the optimization process. The numerical data is presented in dimensionless magnitudes, where the horizontal axis represents the values f / f 1 being f 1 is the first natural frequency, while the vertical axis represents the values FRF/FRF( f 1 ) being FRF( f 1 ) the value of the FRF at each response point for the first natural frequency. The results show that the optimization process led to a good fitting between the numerical and experimental results.
This research deals with the disagreements that can be found between the dynamic response from FE simulations of assembled aircraft structures and experimental testing, usually carried out through a Ground Vibration Test (GVT). These disagreements can be noticed when comparing the Frequency Response Functions (FRF) obtained from both methods, and they are usually due to the complex behaviour of the connections involved, which are difficult to simulate in the FE model. The mechanical properties of these connections are included in the global stiffness matrix of the assembly. Hence, this work presents a novel formulation of the dynamic equilibrium equation using modal superposition with a truncated and constant modal base. This formulation isolates such joint mechanical properties from the rest of the assembly in order to focus on them. The accuracy of this formulation has been tested
119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE
AID:4121 /FLA
[m5G; v1.221; Prn:27/07/2017; 8:46] P.12 (1-12)
S. Hernández et al. / Aerospace Science and Technology ••• (••••) •••–•••
12
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
Fig. 26. Fitting of the experimental results and the numerical FRF at three RP in parallel. Response in RP3.
18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44
in a FE model that aims to represent the connection between the rear part of a fuselage and the tail cone of an aircraft, providing satisfactory results. Furthermore, an optimization approach has been developed based on the aforementioned formulation in order to obtain the mechanical properties of the joints that obtain the best correlation between the FRF of the FE model with the targeted results. The objective function tries to minimize the differences between multiple target and numerical response points obtaining adjustments of good quality. Once the optimized parameters are obtained they can be incorporated to the numerical model to recalculate the modal base and the eigenvalues, and this process can be done iteratively until a good convergence is achieved. The results obtained from this optimization approach show relevant improvements in the correlation of both FRFs. It can be highlighted that as the dynamic response of the structures is tested in several response points, a large number of FRFs may be analysed. When the number of FRF increases, the response of all the FRF cannot be fitted completely. However, since the objective function of the optimization approach considers the absolute error between the FRFs the optimization algorithm will focus on the response points where the discrepancies are higher, leading to an efficient global improvement in the dynamic behaviour the FE model. The methodology has been applied to a FE model that aims to represent a connection between two fuselage sections, and has been validated using a FE model of a real aircraft system for which experimental data is available.
45 46
Conflict of interest statement
47 48
None declared.
49 50
Acknowledgements
51 52 53 54 55 56
The research leading to these results is part of a collaboration research project between Airbus Operations and University of Coruña. The authors also acknowledge funding received from the Galician Government through research grant GRC2013-056 (including FEDER funding).
57 58
References
59 60 61 62 63 64 65 66
[1] W.C. Hurty, Vibrations of structural systems by component mode synthesis, J. Eng. Mech. 86 (1960) 51–69. [2] R.R. Craig, M.C.C. Bampton, Coupling of substructures for dynamic analysis, AIAA J. 6 (1968) 1313–1319. [3] D.-M. Tran, Component mode synthesis methods using interface modes. Application to structures with cyclic symmetry, Comput. Struct. 79 (2001) 209–222. [4] K.J. Bathe, J. Dong, Component mode synthesis with subspace iterations for controlled accuracy of frequency and mode shape solutions, Comput. Struct. 139 (2014) 28–32.
[5] B. Jetmundsen, R.L. Bielawa, W.G. Fannelly, Generalized frequency domain substructure synthesis, J. Am. Helicopter Soc. 33 (1988) 54–64. [6] W. Liu, Structural Dynamic Analysis and Testing of Coupled Structures, PhD thesis, Department of Mechanical Engineering, Imperial College London, 2000. [7] W. Liu, D.J. Ewins, Substructure synthesis via elastic media, J. Sound Vib. 257 (2) (2002) 361–379. [8] A. Sestieri, Structural dynamic modification, Sadhana. Acad. Proc. Eng. Sci. 135 (3) (2000) 247–259. [9] E. Menga, Structural joint identification methodology, in: Structures under Shock and Impact XII, 2012. [10] E. Menga, S. Hernández, S. Moledo, C. López, Nonlinear dynamic analysis of assembled structures, in: 16th International Conference on Aeroelasticity and Structural Dynamics, IFASD2015, St. Petersburg, Russia, 2015. [11] R.A. Ibrahim, C.L. Pettit, Uncertainties and dynamic problems of bolted joints and other fasteners, J. Sound Vib. 279 (2005) 857–936. [12] S. Noll, J.T. Dreyer, R. Singh, Identification of dynamic stiffness matrices of elastomeric joints using direct and inverse methods, Mech. Syst. Signal Process. 39 (2013) 227–244. [13] A.B. Stanbridge, K.Y. Sanliturk, D.J. Ewins, J.V. Ferreira, Experimental investigation of dry friction damping and cubic stiffness nonlinearity, in: ASME Design Technical Conference, Pittsburgh, Pennsylvania, 2001. [14] H. Jalali, H. Ahmadian, J.E. Mottershead, Identification of nonlinear bolted lap-joints parameters by force-state mapping, Int. J. Solids Struct. 44 (2007) 8087–8105. [15] Ö. Arslan, M. Aykan, H.N. Özgüven, Parametric identification of structural nonlinearities from measured frequency response data, Mech. Syst. Signal Process. 25 (2011) 1112–1125. [16] H. Ahmadian, H. Jalali, Identification of bolted lap-joints parameters in assembled structures, Mech. Syst. Signal Process. 21 (2007) 1041–1050. [17] M. Iranzad, H. Ahmadian, Identification of nonlinear bolted lap joint models, Comput. Struct. 96–97 (2012) 1–8. [18] J. Kim, J.C. Yoon, B.S. Kang, Finite element analysis and modelling of structure with bolted joints, Appl. Math. Model. 31 (2007) 895–911. [19] D.J. Segalman, Modelling joint friction in structural dynamics, Struct. Control Health Monit. 13 (2006) 430–453. [20] W.P. Petrov, D.J. Ewins, Analytical formulation on friction interface elements for analysis of nonlinear multi-harmonic vibrations of bladed discs, J. Turbomach. 125 (2003) 364–371. [21] E.P. Petrov, A method for use of cyclic symmetry properties in analysis of nonlinear multiharmonic vibrations of bladed disks, J. Turbomach. 126 (2004) 175–183. [22] H. Ahmadian, H. Jalali, Generic element formulation for modelling bolted lap joints, Mech. Syst. Signal Process. 21 (2007) 2318–2334. [23] F. Wei, G.T. Zheng, Nonlinear vibration analysis of spacecraft with local nonlinearity, Mech. Syst. Signal Process. 24 (2010) 481–490. [24] S. Huang, E.P. Petrov, D.J. Ewins, Comprehensive analysis of periodic regimes of forced vibration for structures with nonlinear snap-through springs, Appl. Mech. Mater. 5–6 (2006) 3–13. [25] M. Oldfield, H. Ouyang, J.E. Mottershead, Simplified models of bolted joints under harmonic loading, Comput. Struct. 84 (2005) 25–33. [26] T. Kalaycioglu, H.N. Özgüven, Nonlinear structural modification and nonlinear coupling, Mech. Syst. Signal Process. 46 (2014) 289–306. [27] D. Göge, Automatic updating of large aircraft models using experimental data from ground vibration testing, Aerosp. Sci. Technol. 7 (2003) 33–45. [28] J.E. Mottershead, M. Link, M.I. Friswell, The sensitivity method in finite element model updating: a tutorial, Mech. Syst. Signal Process. 25 (2011) 2275–2296. [29] S.V. Modak, T.K. Kundra, B.C. Nakra, Comparative study of model updating methods using simulated experimental data, Comput. Struct. 80 (2002) 437–447. [30] B. Peeters, P. De Baets, L. Mosenich, A. Vecchio, H. Van der Auweraer, F. Lambert, Ground vibration testing in the aeroelastic design and certification of a small composite aircraft, in: 46th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics & Materials Conference, Austin (TX), USA, 2005. [31] G. Groll, Windmilling in Aero-Engines, PhD thesis, Department of Mechanical Engineering, Imperial College London, 2000. [32] G. Gloth, M. Sinapius, Analysis of swept-sine runs during modal identification, Mech. Syst. Signal Process. 18 (2004) 1421–1441. [33] S. Hernández, E. Menga, A. Baldomir, C. López, M. Cid, S. Moledo, A methodology for identification of dynamic parameters in assembled aircraft structures, WIT Trans. Model. Simul. 55 (2013) 353–365. [34] K. Shingu, Y. Aoki, T. Irie, K. Mitsui, K. Ogawa, F. Nanaumi, A study on damping characteristics of shell and spatial structures – damping ratios of a cylindrical shell, in: 13th World Conference on Earthquake Engineering, Vancouver, Canada, 2004. [35] K. Shingu, K. Hiratsuka, K. Mitsui, T. Kawashima, N. Kondo, K. Ogawa, T. Irie, M. Yukawa, T. Otsuka, Study on vibration damping characteristics of a spherical shell and damping tendency of shell and spatial structures, in: 15th World Conference on Earthquake Engineering, Lisbon, Portugal, 2013. [36] MSC NASTRAN Documentation. [37] MATLAB R2012b Documentation.
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
[m5G; v1.221; Prn:27/07/2017; 8:46] P.1 (1-12)
Aerospace Science and Technology ••• (••••) •••–•••
1
67
Contents lists available at ScienceDirect
2
68
3
69
Aerospace Science and Technology
4
70
5
71
6
72
www.elsevier.com/locate/aescte
7
73
8
74
9
75
10
76
11 12 13
Optimization approach for identification of dynamic parameters of localized joints of aircraft assembled structures
14 15
77 78 79 80
a,∗
b
a
a
a
a
81
17
S. Hernández , E. Menga , S. Moledo , L.E. Romera , A. Baldomir , C. López , M. Cid Montoya a
18
a
84
16
19
Structural Mechanics Group, School of Civil Engineering, Universidade da Coruña, A Coruña, Spain b Component Loads and Aeroelastic Department – EGLRG, Airbus Operations Getafe, Getafe, Spain
20
25 26 27
i n f o
a b s t r a c t
88 89
Article history: Received 12 September 2016 Received in revised form 25 May 2017 Accepted 17 July 2017 Available online xxxx
28 29 30 31 32 33
85 87
a r t i c l e
23 24
83
86
21 22
82
Keywords: Aircraft design Assembled structures Ground vibration test Dynamic analysis Optimization methods
34 35 36 37 38 39 40 41 42
Ground Vibration Test (GVT) is a crucial step in the design process of a new aircraft. The test involves the analysis of different load configurations that helps to understand the real structural behaviour, to certificate some extreme load cases and also to validate and improve the Finite Element (FE) model used in the design. The validation of the FE model is done by making a comparison between Frequency Response Functions (FRF) from the GVT and the FE model at some response points. Ideally, numerical and experimental results should be in agreement. This is certainly true in the case of single components but in the case of assemblies it is common that numerical and experimental results show discrepancies, especially in the case of dynamic analysis. These differences increase as the number of components in the assembly grows, because connections are usually not represented sufficiently well in the FE model. The components of these assembled structures are usually connected at a small number of points, and they can be seen as local sources of stiffness and damping. In the FE model these joints can be represented as lumped spring/damper elements in a simplified way. In this paper a formulation of the dynamic equilibrium equation employing a decoupled approach for the joint mechanical properties of an assembly is presented. Later, an optimization process based on this formulation that involves the joint properties as design variables has been developed in order to identify the joint parameters that improve the adjustment between numerical and experimental FRF in the case of harmonic excitations. The approach has been applied to a FE model that represents the connection between a segment of fuselage and the tail cone and to a real aircraft system where experimental data is available. The results achieved with this procedure are very efficient and promising. © 2017 Elsevier Masson SAS. All rights reserved.
90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108
43
109
44
110
45 46
111
1. Introduction
47 48 49 50 51 52 53 54 55 56 57 58 59
In the manufacturing process of an aircraft its components are built individually, even at different locations. Wings, fuselage sections, vertical or horizontal stabilizers and other components are assembled after their production in a final assembly line. These components are usually connected at few locations to create an assembly. Fig. 1a shows an example of a single component of an aircraft fuselage while Fig. 1b presents an assembly of two components. In the case of single components, modelling and dynamic prediction has been well developed, and accurate predictions may be achieved. However, in the case of assemblies the quality of the results decreases with the complexity of the connection. This arises
60 61 62 63 64 65 66
*
Corresponding author. E-mail address: [email protected] (S. Hernández).
http://dx.doi.org/10.1016/j.ast.2017.07.026 1270-9638/© 2017 Elsevier Masson SAS. All rights reserved.
due to two possibilities: the component model is not sufficiently accurate or the connection has not been represented satisfactorily. The second reason is frequently the most probable. The appropriate characterization of mechanical properties in assemblies is a recurrent problem for both static and dynamic analyses. Fitting the correct value of stiffness and local damping involved in this type of connections is crucial to obtain accurate models able to predict reliable responses to any kind of loads. To efficiently analyse the dynamic behaviour of assembled structures several methods using substructure coupling have been developed. The Component Mode System (CMS) method was proposed by Hurty [1], and has been widely used and developed [2–4] for the last forty years to predict the frequencies and modal behaviour of complex FE models on the basis of the subsystems modal properties that take part of them. In the case of using experimental Frequency Response Functions (FRF) as response data, the FRF based coupling method proposed by Jetmundsen [5] en-
112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE
2
AID:4121 /FLA
[m5G; v1.221; Prn:27/07/2017; 8:46] P.2 (1-12)
S. Hernández et al. / Aerospace Science and Technology ••• (••••) •••–•••
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39
Fig. 1. Example of components and assembled structures.
40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
ables the structural synthesis in the frequency domain. Both methods allow combining experimental and numerical results and were initially developed considering rigid joints connections between assembled parts. In order to improve the dynamic analysis of assembled structures considering more realistic joint connections, several modifications of these two methods have been proposed, like the General Joint Description Method (GJDM), to take into account the effect of elastic media between subsystems [6,7]. Other approach to study the dynamic response of assembled structures is the Structural Dynamic Modification Method (SDMM), which allows evaluating the incidence of a set of modifications in the dynamic response of the structure without the need to perform FE simulations repeatedly. The SDMM is particularly efficient when the modifications are concentrated in lumped elements. The modified structural response can be expressed as a function of the baseline FE model dynamic database and the set of modifications. A good theoretical approach of SDMM and some work based on it is presented in [8–10]. From a modelling point of view connections can be seen as local sources of stiffness and damping, in most cases with a nonlinear response. Connections are usually performed with prestressed bolts, rivets, pins, adhesives or elastomeric parts, and the complex behaviour of these elements with contact forces, friction, precompression loads and nonlinear material laws substantially modifies the dynamic response of the global system due to changes in stiffness and damping [11–13]. The accurate characterization and mod-
elling of the dynamic behaviour of these connections is a prolific field of research and important efforts have been made in both the numerical modelling and parameter identification by comparing with experimental results [14–17]. They are usually defined by nonlinear constitutive laws, which parameters might depend from a great variety of aspects, such as load application points, load frequency or magnitude, relative displacement value or load duration. Detailed three-dimensional joint models are accurate but not applicable in the analysis of large structures with a high number of connections [18]. The most commonly used models are lumped parameters with the joint modelled as a single point [19–21] and to a lesser extent thin layer elements [22]. Obtaining nonlinear dynamic responses is a difficult task, especially for large structures, due to the necessity of employing iterative methods in time or frequency domain [23] and the complexity of nonlinear vibrations. Among the solutions currently employed to solve this problem, some can be highlighted such as the Harmonic Balance Method (HBM) in the frequency domain [22,24]. HBM is a general method based on identifying periodic solutions as truncated Fourier series with a variable number of harmonics; simplified models for nonlinear time history responses [25], or nonlinear internal force vectors expressed as an equivalent stiffness matrix with dependency of vibration amplitude [26]. During the project of a new aircraft a Finite Element (FE) model is created based on expert knowledge, FE modelling state-of-theart and information learnt from previous successful designs. Once all parts of the new aircraft are designed and manufactured a Ground Vibration Test (GVT) is performed. The GVT is done for validating and improving the structural model. During the GVT several shakers or Driving Points (DP) are located at specific positions of the structure in order to excite the full aircraft, and a great amount of measurement devices, like accelerometers or Response Points (RP) are located throughout the whole aircraft to register as much information as possible. Fig. 2 shows a GVT carried out in a civil aircraft and in a military aircraft. The GVT is carried out using different load configurations. Meanwhile the GVT is performed the FE model is updated. What is learnt from the experiment is introduced in the model, and predictions from the model are used to obtain better experimental results. Since the test is performed in the final stages after the assembly of the prototype, one of the priorities is to perform testing and model updating in the shortest period of time. Updating the mathematical model is crucial to make reliable flutter predictions for flight tests [27–29], and this is why recognizing experimentally the low-frequency modes of the whole aircraft structure is the first scope [30]. It is also an objective of GVT the certification of Fan-Blade-Off and Windmilling analysis [31], because these extreme cases can never be tested under real conditions. Additional testing is required in the case of complex aircraft designs trying to quantify the non-linear behaviour of the structure comparing the obtained responses with different load levels. In the last years new issues have arisen like passenger’s comfort. In the past the use of the Phase Resonance Method or socalled Normal Mode Testing has been used for GVT of large aircrafts, which essentially consist of applying single sine excitations at the structural natural frequencies. This method is quite good for separating closely spaced modes, but has a very time-consuming testing procedure. Hence, it has been substituted by the so-called phase separation techniques that find the aircraft modes by evaluating the FRFs. The signals used in this method are harmonic, periodic signals like multi-sines, swept-sines [32] or periodic chirp, transient signals like impulses or random. The main target of matching FE models and GVT is to validate the FE model in order to predict the behaviour of the structure. However it is not uncommon to find many discrepancies between both results. For many years the lack of reliability when modelling
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE AID:4121 /FLA
[m5G; v1.221; Prn:27/07/2017; 8:46] P.3 (1-12)
S. Hernández et al. / Aerospace Science and Technology ••• (••••) •••–•••
3
where m, c, and k are the mass, modal damping and stiffness matrix, while u and f represent the displacement and the external forces vectors. In the case of assembled structures, the stiffness matrix represents the stiffness of the whole assembly. Since this work pursues to improve the fitting of the FE results against experimental testing by obtaining the adequate values of the mechanical properties of the joints that connect the components of the assembly, it would be necessary to isolate such properties from the dynamic equilibrium equation. In this sense, a novel formulation with the aim of identifying the dynamic parameters of the joints in assembled structures in order to fit numerical FRF with experimental FRF is developed in this section for being applied later on. As seen in Hernández [33], assuming linear spring joints, the stiffness matrix k can be divided in two parts (ks + kns ) based on the SDMM:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 18 19 20 21 22 24 25 26 27 28 29 31
(2)
where the matrix k g has the products of the structural damping by the variable spring stiffness:
32 33
kg =
35
⎛
Fig. 2. Ground Vibration Test of a civil and a military aircraft.
37
41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59
the junctions in complex assemblies has been underestimated. This research presents a decoupled formulation of the dynamic equilibrium equation that isolates the local stiffness and damping of the FE model connections from the rest of the structure, aiming to identify through an optimization process their optimum values in order to improve the performance of FE models of assembled structures for reaching a good agreement between experimental and numerical results. The novelty of this research lies in the formulation exposed in Section 2, which was developed in the framework of an industrial collaboration with Airbus Operations for improving the correlation between the experimental and numerical results of assembled structures connected at a small number of joints. In this sense, the dynamic behaviour of the parts of the assembly is assumed to be well represented in the FE analyses, being the focus of the discrepancies in the numerical behaviour of the joints, which is certainly the most usual case in real situations. This methodology is successfully applied to a simplified FE model that represents the assembly between a rear fuselage section and the tail cone of an aircraft and to an aircraft system that consist of a Vapour Cycle Refrigeration Unit (VCRU) connected to the fuselage of an aircraft. The latter case is studied considering real experimental data.
60 61 62 63 64
2. Formulation of the dynamic analysis of assembled structures through a decoupled approach The traditional formulation of the dynamic equilibrium of a structure is:
65 66
71 72 73 74 75 76 77 78 79 80 81 83 84 85 86 87 89 90 91 92 93 94 96 97 98 99
34
40
70
95
k = (kns + ks )(1 + ig ) = kns + ks + ik g
30
39
69
88
In a limit case, the entire spring stiffness might be assigned to ks , while kns would contain just the stiffness of the separate components. However, in practical applications this would worsen the behaviour of the assembly since the components will not even be partially connected. Considering structural damping g, the total structural stiffness can be split into the sum of three matrices.
23
38
68
82
– A matrix kns that contains the stiffness of all the separate parts that builds the assembly and a partial and constant part of the spring stiffness. – A matrix ks than contains the rest of the uncoupled spring stiffness, considering it as variable.
17
36
67
¨ + cu˙ + ku = f mu
(1)
100
.. ⎜ . ⎜. . . (k + k ) g nsi si i ⎜ ⎜ −(k + k ) g nsi si i ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ . . .
101
⎞
. . .
⎟ ⎟ ... p −(knsi + ksi ) g i . . . ⎟ ⎟ ...q (knsi + ksi ) g i ⎟ ⎟ . . ⎟ . . .. . ⎟ . . ⎟ . . . (knsj + ksj ) g j −(knsj + ksj ) g j . . . ⎟ ⎟ ...r ⎟ ...s −(knsj + ksj ) g j (knsj + ksj ) g j ⎟ ⎠ . .. . . . . . . . . . . . .
p
q
r
s
(3) In which p and q, r and s, are the degrees of freedom that connect the springs i and j respectively, as seen in Equation (3). The total stiffness of the joint, represented as ki is the sum of k si and knsi , and g i is the structural damping value of the joint, which is the local damping parameter of the joint. The matrix ks has the variable stiffness information of the joints:
⎛.
⎜. . . k −ksi ⎜ si ⎜ −ksi k si ⎜ ⎜ .. ⎜ ks = ⎜ . ⎜ ⎜ ⎜ ⎝ .. .
.. .
p
q
... .. . . −ksj . . . ksj −ksj k sj .. . .. .. . . ..
r
s
⎟ ... p ⎟ ⎟ ...q ⎟ ⎟ ⎟ ⎟ ⎟ . . .⎟ . . .r ⎟ ⎠ ...s .. .
103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
⎞
.. .
..
102
121 122 123 124
(4)
125 126 127 128 129 130 131 132
JID:AESCTE
AID:4121 /FLA
[m5G; v1.221; Prn:27/07/2017; 8:46] P.4 (1-12)
S. Hernández et al. / Aerospace Science and Technology ••• (••••) •••–•••
4
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
where p and q are the joints that define the i spring. Thus, the dynamic equilibrium under a vector of forces f can be rewritten as:
20 21 22
Defining the displacement vector u as a function of the spectral matrix which contains the natural vibration modes, u becomes:
u = q,
(6)
where is the spectral matrix that contains the eigenvectors associated with the matrix kns normalized to the total mass. The vectors f and q can be written as
q = ve i ωt .
(7)
Introducing these expressions in (5) the dynamic equilibrium becomes
Iq¨ + Cq˙ + Kns q + Ks q + iK g q = F,
23
C = Φ cΦ
24
and
25 26 27 28 29
32
T
F = Φ p,
Kns = Φ T kns Φ
⎛
35 36 37 38 39 40 41 42 43 44 45 46 47 48
Ks = Φ T ks Φ
⎜
ω12 . . .
. Kns = ⎝ ..
..
K g = Φ T k g Φ,
⎞
⎛
.. ⎟ . ⎠
C=⎝
⎜
. 2 . . . ωm
33 34
(9)
(10)
where the matrices C and Kns are easily built from the modal information as follows,
30 31
(8)
where T
2ζ1 ω1
.. .
⎞
... ⎟ .. .. ⎠ , (11) . . . . . 2ζm ωm
being m the number of natural frequencies. The vectors p, F and v in general are complex, so
p = pr + ipi
F = Fr + iFi
v = vr + ivi .
Therefore, it can be written that
Rns + Ks ωC + K g
−(ωC + K g ) Rns + Ks
vr vi
=
Fr Fi
(12)
(13)
,
where
Rns = Kns − ω2 I.
(14)
Solving the system of linear equations (13) the vector of dis¨ can be obtained placements u and the vector of accelerations u for each type of dynamic load as
49 50 51 52 53 54 55 56 57 58 59 60 61
u = q = |vr + ivi |e
64 65 66
¨ = −ω |vr + ivi |e u 2
i ωt
.
(15)
A specific displacement uk or acceleration u¨ k will be written as
uk = ψ k vr + i (ψ vi ) e i ωt
u¨ k = −ω ψ k vr + i (ψ vi ) e i ωt , (16) 2
where ψ k is the k-th row of matrix . This can be expressed as,
uk = (ψ k vr )2 + (ψ k vi )2 e i (ω+θk )t = Ae i (ω+θk )t u¨ k = −ω2 (ψ k vr )2 + (ψ k vi )2 e i (ω+θk )t = −ω2 Ae i (ω+θk )t
(17) (18)
where θk is the phase angle with a value of
62 63
i ωt
θk = arctg
ψ k vi ψ k vr
69
(5)
18 19
68 70
¨ + cu˙ + ku = mu¨ + cu˙ + (kns + ks + iks )u = f. mu
f = pe i ωt
67
(19)
and therefore A is the magnitude of the FRF of displacement uk for the load with frequency ω .
71
Fig. 3. Springs definition.
This formulation allows to obtain the displacements, velocities and accelerations of a structure in a novel way. It is only required the spectral matrix ns , without the variable springs stiffness, and the natural vibration frequencies. The main advantage is its flexibility. From the global modal information, which only needs to be obtained once, several fast trials with different springs definition can be made by just changing the matrices k g and ks . The analysis of one FRF can be done following these steps: – Analyse the eigenvalue problem and build the spectral matrix ns , after removing the variable spring stiffness. – Assemble of stiffness matrices k g and ks . – Compute K g and Ks as seen in expression (10). – Create the matrix Rns as seen in expression (14). – Solve the system of equations in (13), and compute expression (16) to obtain the displacements and accelerations. The final purpose of this formulation is to produce a fitting between numerical and experimental results. The main advantage is the fact that the truncated spectral matrix ns remains constant, which simplifies the process, removing the need of solving the dynamic modal analysis of the full structure every time a modification in the joints is made. These modifications affect only the matrices k g and ks .
72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98
3. Validation of the uncoupled formulation
99
The formulation proposed to solve this problem divides the stiffness matrix into the sum of three as seen in expression (2). In order to validate this formulation several numerical tests with different ratio of uncoupled springs stiffness have been made. The analysis has been done in the structural laminar model defined in Fig. 3, which represents the connection between two sections of an aircraft. The model consists of an 8 meters long cylinder of 4 meters of radius and a truncated cone of 4 meters long and varying radius from 4 to 2 meters. The cylinder that represents the fuselage is considered built-in and a sinusoidal unitary load is applied at the tip of the cone that represents the tail of the fuselage, as shown in Fig. 4. A constant value of 2 mm has been applied to the shell elements of the model and aluminium, with Young’s modulus value of 70000 MPa and Poisson’s ratio of 0.2, has been selected. The connection between the cylinder and the cone is composed by 5 joints, with stiffness in 9 degree of freedom in total. The position of the joints can be identified by the angle α as presented in Fig. 5, which takes the following values: 45◦ , 135◦ , 225◦ , 270◦ and 315◦ . The location and connectivity of the joints can be seen as well in Fig. 5. These joints represent the assembly between a segment of fuselage and the tail cone of a real aircraft. In order to validate the formulation several analyses with different coupled stiffness ratio have been carried out for coupled stiffness in range 40% to 100%. For instance, if a ratio of 60% of coupled stiffness is defined it means that 40% of the stiffness of the lumped elements is in the matrix ks (matrix of variable springs stiffness). The remaining stiffness is obtained from the modal analysis of the finite element model changing the stiffness of the joint spring elements to the 60% of their initial value. This stiffness is computed directly using the matrix Kns with no more information than the eigenvalues as seen in (11).
101
100 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE AID:4121 /FLA
[m5G; v1.221; Prn:27/07/2017; 8:46] P.5 (1-12)
S. Hernández et al. / Aerospace Science and Technology ••• (••••) •••–•••
5
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11 12
77
Fig. 4. Model geometry and load location.
78
13
79
14
80
15
81
16
82
17
83
18
Fig. 7. Detail of the comparison of the uncoupled formulation at node 5432 DOF 2, in the frequency range from 15 to 19 Hz, without considering structural damping.
19 20
The formulation shows accurate results at different ratios of coupled stiffness when compared with the 100% coupled solution obtained in MATLAB. Moreover, a validation of the formulation was accomplished after comparing the results of the 100% coupled solution with the numerical one obtained with MSC NASTRAN through a SOL111, a modal frequency response analysis. In Fig. 7, which is a detail of Fig. 6 in the frequency range from 15 to 19 Hz, the effects of modifying the ratio of coupled stiffness can be recognized. It is clear that with higher values of coupled stiffness better results will be reached, being the 100% case of coupled stiffness the exact solution. Fig. 8 and Fig. 9 are the same as Fig. 6 and Fig. 7 but considering added structural damping with a value of 0.1. Once again the comparison with MSC NASTRAN has been done, and a total superposition was obtained. Even at low levels of coupled stiffness the formulation shows to be quite precise.
22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
Fig. 5. Location and connectivity of joints.
38 39 40 41 42 43 44 45 46 47 48 49 50
85 86
21
37
84
The structure has a constant modal damping of 2%, whereas each spring has a stiffness value of 5000 N/m. The damping ratio of 2% is based on experimental vibration studies found in literature [34,35]. Modal analysis has been performed with MSC NASTRAN [36] to obtain the spectral matrix and the natural frequencies on a range from 0 to 100 Hz, which included 229 modes. The formulation developed in the previous section has been implemented in MATLAB [37]. The results shown in Figs. 6 to 9 were obtained with a vertical and unitary harmonic force located at node 5432 (Fig. 4) and the response studied was the vertical acceleration at the same node.
87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103
4. Identification of mechanical parameters by an optimization process
104
In assembled structures certain parameters such as the stiffness or the local damping related to joints may be unknown. In this research the fitting between experimental results obtained from the GVT and numerical results from the FE model has been done by identifying these parameters. Thanks to the uncoupled formulation any change can be easily made in the stiffness matrix without any further information from the FE model. Just varying the uncoupled matrix ks a better fitting between experimental information and numerical results can be achieved.
107
105 106 108 109 110 111 112 113 114 115 116
51
117
52
118
53
119
54
120
55
121
56
122
57
123
58
124
59
125
60
126
61
127
62
128
63
129
64
130
65 66
131
Fig. 6. Comparison of the uncoupled formulation at node 5432 DOF 2 without considering structural damping.
132
JID:AESCTE
AID:4121 /FLA
[m5G; v1.221; Prn:27/07/2017; 8:46] P.6 (1-12)
S. Hernández et al. / Aerospace Science and Technology ••• (••••) •••–•••
6
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
Fig. 8. Comparison of the uncoupled formulation at node 5432 DOF 2 considering a structural damping of 0.1.
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27 29 30 31 32 33 34 35
Fig. 9. Detail of the uncoupled formulation at node 5432 DOF 2, in the frequency range from 15 to 19 Hz, considering a structural damping of 0.1.
36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
93
Fig. 10. Strategy of the optimization process.
28
Furthermore, an optimization process can be defined to find the best value of the objective function. The design variables of the problem are the parameters involved in the connection (spring stiffness and structural damping) and also a set of modal damping values that cover the ranges of the frequencies involved. The objective function for this optimization problem was defined as the sum of two terms, one related to the natural frequencies and another related to the FRF value at those frequencies. Minimization of differences between the numerical and experimental natural frequencies is done in the first term of the objective function, approaching horizontally the numerical and experimental FRF. The second term of the objective function aims to fit the FRF vertically by minimizing the difference between the peaks values, which are the FRFs evaluated at their natural frequencies. The objective function (20) is formulated as the sum of squares of the differences between numerical i-th natural frequency ωiN U and experimental i-th natural frequency ωi E X plus the sum of squares of the differences between numerical values of the FRF evaluated at the i-th natural frequency, FRF N U (ωiN U ), and the experimental value of the FRF evaluated at the i-th natural frequency, FRF E X (ωi E X ).
is the i-th experimental natural frequency, and ωiN U is the i-th numerical natural frequency. In expression (20) the experimental values extracted from the GVT do not vary since they are the target values of the optimization process. The remaining terms are evaluated during the optimization process, because any change in a parameter would change the numerical FRF, not only the position of the natural frequencies but also the peak values at those frequencies. Fig. 10 shows how the fitting between numerical and experimental is reached. The numerical natural frequencies are defined as the points where the maximum values of the FRFs are located. To start the optimization process these natural frequencies are obtained from the solution of the eigenvalue problem. Once all numerical natural frequencies are known a range of frequencies ω is defined in order to reduce the analysis time and evaluate the FRF only around the desired frequencies. Initially, ranges are defined considering the numerical and experimental natural frequencies. The optimization process assures that any changes in the design variables will produce numerical results that are included within the interval [min( w iN U − w , w i E X − w ), max( w iN U + w , w i E X + w )]. Expression (20) is formulated for a single FRF. It might be more suitable to fit simultaneously as many FRFs of degrees of freedoms as possible in order to enhance the accuracy between numerical and experimental results. A more general objective function that includes more response points, and consequently more FRFs, can be defined as
F=
n
F= (ωi E X − ωiN U )2
N FRF
+
2 FRF E X (ωi E X ) − FRF N U (ωiN U )
N w
λj (ωi j E X − ωi jN U )2
j =1
i =1 n
+ (20)
i =1
where n is the number of peaks considered in the optimization, FRF E X is the experimental FRF, FRF N U is the numerical FRF, ωi E X
94
i =1
Nw
95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125
2 FRF j E X (ωi j E X ) − FRF jN U (ωi jN U ) ,
126
(21)
i =1
where N FRF is the number of FRF taken into account, λ j is a weight factor that has been set to 1 in this study, N w is the number of frequencies considered, FRF j E X is the j-th experimental FRF, FRF jN U
127 128 129 130 131 132
JID:AESCTE AID:4121 /FLA
[m5G; v1.221; Prn:27/07/2017; 8:46] P.7 (1-12)
S. Hernández et al. / Aerospace Science and Technology ••• (••••) •••–•••
7
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14 15
80
Fig. 11. Fitting of FRF of vertical acceleration at node 5432.
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31 32
97
Fig. 12. Location of the parameters taken into account in the optimization case with 4 FRF response points simultaneously.
33 34 35 36
99
is the j-th numerical FRF, ωi j E X is the i-th experimental natural frequency of the j-th FRF, and ωi jN U is the i-th numerical natural frequency of the j-th FRF.
37 38 39
5. Application example modelling an assembled aircraft structure
40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The formulation exposed previously has been applied to the structural model defined in Fig. 4, which represents the connection of two fuselage segments of an aircraft. Fig. 5 shows the connectivity of this section, composed by 5 springs working in 9 DOF in total. The experimental data has been simulated by changing some parameters of the model, such as the springs stiffness or the modal damping. Then a frequency domain analysis has been carried out using MSC NASTRAN. The results of this analysis are used as the objective ones in the optimization process, simulating the experimental results obtained from a GVT. Two different configurations have been taken into account in this model: one considering as design variables the spring’s stiffness and the modal damping of several natural frequencies and another considering in addition as design variables the local damping modelled as structural damping. The optimization problem was solved using MATLAB with a Sequential Quadratic Programming (SQP) algorithm. In both cases the dynamic excitation was a harmonic unitary force applied in vertical direction at node 5432, on the top tip of the cone (Fig. 4).
61 62 63
5.1. Optimization without including structural damping as design variable
64 65 66
98
Firstly, an optimization problem with only one response point was solved. This problem uses the objective function defined in
expression (20). The optimization process was formulated considering as design variables the stiffness of the 5 springs, which are 9 design variables (one for each DOF), plus 9 modal damping coefficients selected at the natural frequencies excited in this response. The response considered was the FRF of the vertical acceleration at node 5432. The spring stiffness values of the initial design are a 90% of the spring stiffness values employed in the objective design. Moreover, the values of the damping were modified from a 2% in the objective design to a 3% in the initial design. The objective design in this academic example is obtained numerically through NASTRAN, however in practical applications it comes from experimental tests. These values of initial and objective parameters were used in all cases studied in this research. The number of eigenvectors used is the same as in the validation example. The proposed optimization process is able to reduce the objective function from an initial value of 51.8271 to a final value of 0.8668, which means an improvement of 98.32%. The initial, objective and optimized FRF of the response point selected can be seen in Fig. 11. The tolerance established for the objective function was a relative error of 1 × 10−6 in all cases and the algorithm performed 128 iterations before convergence in this example. After seeing the good results achieved by the optimization process a new problem considering 4 different response points was formulated. The final purpose is to use as many response points as possible, which would mean a better accuracy between FE model and experimental results. In this case the selected objective function was the one defined in expression (21). The following RPs were used:
100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128
– Node 5432 DOF 2, which is the vertical direction of the node where the load is applied. – Node 3405 DOF 3, which corresponds to the longitudinal direction of the upper left spring.
129 130 131 132
JID:AESCTE
8
AID:4121 /FLA
[m5G; v1.221; Prn:27/07/2017; 8:46] P.8 (1-12)
S. Hernández et al. / Aerospace Science and Technology ••• (••••) •••–•••
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
Fig. 13. Fitting of vertical acceleration at node 5432 with 4 response points in parallel.
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29 30
95
Fig. 14. Fitting of longitudinal acceleration at node 3405 with 4 response points in parallel.
31 32 33 34 35 36
97
– Node 3443 DOF 1, which corresponds to the horizontal direction of the bottom spring. – Node 3455 DOF 2, which corresponds to the vertical direction of the bottom right spring. The position of the parameters used in this optimization case can be seen in Fig. 12.
37 38 39 40 41 42 43 44 45 46 47 48 49
In the example, whose results are provided from Figs. 13 to 16, 25 design variables have been used; 9 matching the springs stiffness as seen previously in the one parameter optimization case and 16 corresponding to the modal damping coefficients. In this case there are 5 more than in the previous example because some natural frequency might be excited in one parameter and might not be in another one. To achieve good results, all the modal damping coefficients of the modes excited at any response point must be considered. In this problem the results were also very accurate. The initial value of the objective function was 77.1348 and the final value was 5.0214, which means an improvement of 93.49%.
50 51 52
5.2. Optimization including structural damping as design variables
53 54 55 56 57 58 59 60 61 62 63 64 65 66
96
In this example local damping, modelled as structural damping, has been added to the optimization process. Consequently, 9 new design variables were added to the problem, one at each spring component. The initial values of all structural damping are the double of the ones employed in the objective design, which as exposed in section 5.1 represent the experimental data. Following the same idea presented in section 5.1, a simpler optimization considering just one response point was initially formulated. In this case the total number of design variables was 27. The optimization process shows also good results with the structural damping formulation as it can be seen in Fig. 17, in fact the improvement was a 95.05%. On the following pictures g refers to structural damping. The initial point of the optimization process had an ob-
jective function value of 40.2751, and after optimization it was 1.9936. The next step was to carry out the formulation with 4 response points in parallel. This resulted in a bigger optimization problem in which 31 design variables were used. The response points used in this case were the same used in the example of section 5.1. In this example the decrease of the objective function is a bit lower than in previous cases, but as it can be seen in Figs. 18–21 the structural behaviour has drastically been improved. An improvement of 91.97% was achieved, reducing the initial value of the objective function from 59.9372 to 4.8110. It must be pointed out that the response with less magnitude, such as node 3405 and mainly node 3443, have lower improvements that the response points with larger responses (5432 and 3455), maintaining higher discrepancies in comparison with more excited RP. This is due to the use of absolute error in the formulation of the objective function, which was defined intentionally because the main target is to improve the larger responses. However, this can be modified by considering relative instead of absolute errors in the objective function or weighting the term of the objective function associated with the most important responses. In Table 1 there is a summary of the results achieved in this research.
98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121
6. Application to an aircraft system with experimental data available
122 123 124
The methodology presented in Section 4 has been applied to the FE model of an aircraft system for which experimental data is available. The test is carried out to assure an appropriate level of vibrations in the structure. The system is a Vapour Cycle Refrigeration Unit (VCRU), which is usually located in pairs at each side of the fuselage, being connected to it through a suspension system composed by a frame and a set of struts, each one containing a rubber mount that is modelled as a lumped spring
125 126 127 128 129 130 131 132
JID:AESCTE AID:4121 /FLA
[m5G; v1.221; Prn:27/07/2017; 8:46] P.9 (1-12)
S. Hernández et al. / Aerospace Science and Technology ••• (••••) •••–•••
9
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14 15
80
Fig. 15. Fitting of transversal acceleration at node 3443 with 4 response points in parallel.
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31 32
97
Fig. 16. Fitting of vertical acceleration at node 3455 with 4 response points in parallel.
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48 49
114
Fig. 17. Fitting of vertical acceleration at node 5432 with structural damping considering only this response point.
115
50
116
51
117
52
118
53
119
54
120
55
121
56
122
57
123
58
124
59
125
60
126
61
127
62
128
63
129
64
130
65 66
131
Fig. 18. Fitting of vertical acceleration at node 5432 with structural damping and 4 response points in parallel.
132
JID:AESCTE
AID:4121 /FLA
[m5G; v1.221; Prn:27/07/2017; 8:46] P.10 (1-12)
S. Hernández et al. / Aerospace Science and Technology ••• (••••) •••–•••
10
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
Fig. 19. Fitting of longitudinal acceleration at node 3405 with structural damping and 4 response points in parallel.
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
Fig. 20. Fitting of vertical transversal at node 3443 with structural damping and 4 response points in parallel.
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
Fig. 21. Fitting of vertical acceleration at node 3455 with structural damping and 4 response points in parallel.
48 49 50 51
114
Cases
52 53 54
115
Table 1 Summary of all optimizations carried out.
1 RP 4 RP
116
Without structural damping
117
With structural damping
Initial F obj
Final F obj
%
Initial F obj
Final F obj
%
51.8271 77.1348
0.8684 5.0214
98.32 93.49
40.2751 59.9372
1.9936 4.8110
95.05 91.97
55
58 59 60 61 62 63 64 65 66
119 120 121
56 57
118
122
element. The components of the system (compressor, evaporator, fan, . . .) are modelled as concentrated masses that are connected to the frame through rigid elements and springs. Fig. 22 shows the FE model of a VCRU and their location within the fuselage. The response points selected for the experimental validation are the point where the excitation is applied (RP1) and two points of the frame located in the upper side of the frame, near the connection to the fuselage (RP2 and RP3). The locations of the response points are shown in Fig. 23.
The excitation in the experimental test is introduced as a frequency dependent and periodically rotating radial force that is applied at RP1. The accelerations are measured in the three response points at five excitation frequencies. The experimental results provided are a modified version of the results of a real project. In this example, the target is to fit the numerical FRF to the experimental data at the frequency points where information is available. The design variables considered are the modal damping of the structure and the local stiffness and damping of the shock mounts. The modal damping is frequency dependant and changes
123 124 125 126 127 128 129 130 131 132
JID:AESCTE AID:4121 /FLA
[m5G; v1.221; Prn:27/07/2017; 8:46] P.11 (1-12)
S. Hernández et al. / Aerospace Science and Technology ••• (••••) •••–•••
11
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16 17
Fig. 22. FE model of the VCRU and location within the fuselage.
82 83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29 30
95
Fig. 23. Location of the response points selected.
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47 48
Fig. 24. Fitting of the experimental results and the numerical FRF at three RP in parallel. Response in RP1.
Fig. 25. Fitting of the experimental results and the numerical FRF at three RP in parallel. Response in RP2.
49
116
7. Conclusions
51 52
56 57 58 59 60 61 62 63 64 65 66
117 118
53 55
114 115
50
54
113
every 10 Hz in the range of study [0–200] Hz, leading to 21 damping design variables. Figs. 24–26 show the experimental results measured in the three response points at the target excitation frequencies, the initial numerical FRF and the numerical FRF obtained after the optimization process. The numerical data is presented in dimensionless magnitudes, where the horizontal axis represents the values f / f 1 being f 1 is the first natural frequency, while the vertical axis represents the values FRF/FRF( f 1 ) being FRF( f 1 ) the value of the FRF at each response point for the first natural frequency. The results show that the optimization process led to a good fitting between the numerical and experimental results.
This research deals with the disagreements that can be found between the dynamic response from FE simulations of assembled aircraft structures and experimental testing, usually carried out through a Ground Vibration Test (GVT). These disagreements can be noticed when comparing the Frequency Response Functions (FRF) obtained from both methods, and they are usually due to the complex behaviour of the connections involved, which are difficult to simulate in the FE model. The mechanical properties of these connections are included in the global stiffness matrix of the assembly. Hence, this work presents a novel formulation of the dynamic equilibrium equation using modal superposition with a truncated and constant modal base. This formulation isolates such joint mechanical properties from the rest of the assembly in order to focus on them. The accuracy of this formulation has been tested
119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE
AID:4121 /FLA
[m5G; v1.221; Prn:27/07/2017; 8:46] P.12 (1-12)
S. Hernández et al. / Aerospace Science and Technology ••• (••••) •••–•••
12
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
Fig. 26. Fitting of the experimental results and the numerical FRF at three RP in parallel. Response in RP3.
18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44
in a FE model that aims to represent the connection between the rear part of a fuselage and the tail cone of an aircraft, providing satisfactory results. Furthermore, an optimization approach has been developed based on the aforementioned formulation in order to obtain the mechanical properties of the joints that obtain the best correlation between the FRF of the FE model with the targeted results. The objective function tries to minimize the differences between multiple target and numerical response points obtaining adjustments of good quality. Once the optimized parameters are obtained they can be incorporated to the numerical model to recalculate the modal base and the eigenvalues, and this process can be done iteratively until a good convergence is achieved. The results obtained from this optimization approach show relevant improvements in the correlation of both FRFs. It can be highlighted that as the dynamic response of the structures is tested in several response points, a large number of FRFs may be analysed. When the number of FRF increases, the response of all the FRF cannot be fitted completely. However, since the objective function of the optimization approach considers the absolute error between the FRFs the optimization algorithm will focus on the response points where the discrepancies are higher, leading to an efficient global improvement in the dynamic behaviour the FE model. The methodology has been applied to a FE model that aims to represent a connection between two fuselage sections, and has been validated using a FE model of a real aircraft system for which experimental data is available.
45 46
Conflict of interest statement
47 48
None declared.
49 50
Acknowledgements
51 52 53 54 55 56
The research leading to these results is part of a collaboration research project between Airbus Operations and University of Coruña. The authors also acknowledge funding received from the Galician Government through research grant GRC2013-056 (including FEDER funding).
57 58
References
59 60 61 62 63 64 65 66
[1] W.C. Hurty, Vibrations of structural systems by component mode synthesis, J. Eng. Mech. 86 (1960) 51–69. [2] R.R. Craig, M.C.C. Bampton, Coupling of substructures for dynamic analysis, AIAA J. 6 (1968) 1313–1319. [3] D.-M. Tran, Component mode synthesis methods using interface modes. Application to structures with cyclic symmetry, Comput. Struct. 79 (2001) 209–222. [4] K.J. Bathe, J. Dong, Component mode synthesis with subspace iterations for controlled accuracy of frequency and mode shape solutions, Comput. Struct. 139 (2014) 28–32.
[5] B. Jetmundsen, R.L. Bielawa, W.G. Fannelly, Generalized frequency domain substructure synthesis, J. Am. Helicopter Soc. 33 (1988) 54–64. [6] W. Liu, Structural Dynamic Analysis and Testing of Coupled Structures, PhD thesis, Department of Mechanical Engineering, Imperial College London, 2000. [7] W. Liu, D.J. Ewins, Substructure synthesis via elastic media, J. Sound Vib. 257 (2) (2002) 361–379. [8] A. Sestieri, Structural dynamic modification, Sadhana. Acad. Proc. Eng. Sci. 135 (3) (2000) 247–259. [9] E. Menga, Structural joint identification methodology, in: Structures under Shock and Impact XII, 2012. [10] E. Menga, S. Hernández, S. Moledo, C. López, Nonlinear dynamic analysis of assembled structures, in: 16th International Conference on Aeroelasticity and Structural Dynamics, IFASD2015, St. Petersburg, Russia, 2015. [11] R.A. Ibrahim, C.L. Pettit, Uncertainties and dynamic problems of bolted joints and other fasteners, J. Sound Vib. 279 (2005) 857–936. [12] S. Noll, J.T. Dreyer, R. Singh, Identification of dynamic stiffness matrices of elastomeric joints using direct and inverse methods, Mech. Syst. Signal Process. 39 (2013) 227–244. [13] A.B. Stanbridge, K.Y. Sanliturk, D.J. Ewins, J.V. Ferreira, Experimental investigation of dry friction damping and cubic stiffness nonlinearity, in: ASME Design Technical Conference, Pittsburgh, Pennsylvania, 2001. [14] H. Jalali, H. Ahmadian, J.E. Mottershead, Identification of nonlinear bolted lap-joints parameters by force-state mapping, Int. J. Solids Struct. 44 (2007) 8087–8105. [15] Ö. Arslan, M. Aykan, H.N. Özgüven, Parametric identification of structural nonlinearities from measured frequency response data, Mech. Syst. Signal Process. 25 (2011) 1112–1125. [16] H. Ahmadian, H. Jalali, Identification of bolted lap-joints parameters in assembled structures, Mech. Syst. Signal Process. 21 (2007) 1041–1050. [17] M. Iranzad, H. Ahmadian, Identification of nonlinear bolted lap joint models, Comput. Struct. 96–97 (2012) 1–8. [18] J. Kim, J.C. Yoon, B.S. Kang, Finite element analysis and modelling of structure with bolted joints, Appl. Math. Model. 31 (2007) 895–911. [19] D.J. Segalman, Modelling joint friction in structural dynamics, Struct. Control Health Monit. 13 (2006) 430–453. [20] W.P. Petrov, D.J. Ewins, Analytical formulation on friction interface elements for analysis of nonlinear multi-harmonic vibrations of bladed discs, J. Turbomach. 125 (2003) 364–371. [21] E.P. Petrov, A method for use of cyclic symmetry properties in analysis of nonlinear multiharmonic vibrations of bladed disks, J. Turbomach. 126 (2004) 175–183. [22] H. Ahmadian, H. Jalali, Generic element formulation for modelling bolted lap joints, Mech. Syst. Signal Process. 21 (2007) 2318–2334. [23] F. Wei, G.T. Zheng, Nonlinear vibration analysis of spacecraft with local nonlinearity, Mech. Syst. Signal Process. 24 (2010) 481–490. [24] S. Huang, E.P. Petrov, D.J. Ewins, Comprehensive analysis of periodic regimes of forced vibration for structures with nonlinear snap-through springs, Appl. Mech. Mater. 5–6 (2006) 3–13. [25] M. Oldfield, H. Ouyang, J.E. Mottershead, Simplified models of bolted joints under harmonic loading, Comput. Struct. 84 (2005) 25–33. [26] T. Kalaycioglu, H.N. Özgüven, Nonlinear structural modification and nonlinear coupling, Mech. Syst. Signal Process. 46 (2014) 289–306. [27] D. Göge, Automatic updating of large aircraft models using experimental data from ground vibration testing, Aerosp. Sci. Technol. 7 (2003) 33–45. [28] J.E. Mottershead, M. Link, M.I. Friswell, The sensitivity method in finite element model updating: a tutorial, Mech. Syst. Signal Process. 25 (2011) 2275–2296. [29] S.V. Modak, T.K. Kundra, B.C. Nakra, Comparative study of model updating methods using simulated experimental data, Comput. Struct. 80 (2002) 437–447. [30] B. Peeters, P. De Baets, L. Mosenich, A. Vecchio, H. Van der Auweraer, F. Lambert, Ground vibration testing in the aeroelastic design and certification of a small composite aircraft, in: 46th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics & Materials Conference, Austin (TX), USA, 2005. [31] G. Groll, Windmilling in Aero-Engines, PhD thesis, Department of Mechanical Engineering, Imperial College London, 2000. [32] G. Gloth, M. Sinapius, Analysis of swept-sine runs during modal identification, Mech. Syst. Signal Process. 18 (2004) 1421–1441. [33] S. Hernández, E. Menga, A. Baldomir, C. López, M. Cid, S. Moledo, A methodology for identification of dynamic parameters in assembled aircraft structures, WIT Trans. Model. Simul. 55 (2013) 353–365. [34] K. Shingu, Y. Aoki, T. Irie, K. Mitsui, K. Ogawa, F. Nanaumi, A study on damping characteristics of shell and spatial structures – damping ratios of a cylindrical shell, in: 13th World Conference on Earthquake Engineering, Vancouver, Canada, 2004. [35] K. Shingu, K. Hiratsuka, K. Mitsui, T. Kawashima, N. Kondo, K. Ogawa, T. Irie, M. Yukawa, T. Otsuka, Study on vibration damping characteristics of a spherical shell and damping tendency of shell and spatial structures, in: 15th World Conference on Earthquake Engineering, Lisbon, Portugal, 2013. [36] MSC NASTRAN Documentation. [37] MATLAB R2012b Documentation.
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132